DaSEH Project

This was started at the DaSEH Code-A-Thon! Should give some credit to them for coding help.

Analysis Plan from Thesis Proposal

Aim 3 will explore and describe these co-exposures to visually display the county-level spatial distribution of this co-exposure across the United States. Priority areas for public health interventions are areas that have high co-exposures and high vulnerability to environmental hazards like heat. The primary outcome variable for mapping is the co-exposure of heat and glyphosate, defined here as the heat-glyphosate index calculated as 85th percentile high temperature*glyphosate application per square kilometer. I will map this primary outcome on a scale from white to bright red where red is the highest co-exposure. I hypothesize that the Midwest and Great Plains regions will have the highest levels of co-exposure due to the large amounts of glyphosate resistant crops grown there.

Relevant exposures of public health and policy interest include the CDC’s Social Vulnerability Index (SVI), most prevalent type of glyphosate resistant crop, and urban pesticide usage. I will overlay SVI over the glyphosate-heat co-exposure map by outlining counties with high SVI (75th percentile or higher). Geographic unit will be county as this is the unit for glyphosate application rate. I will create separate maps for different crop types (i.e. soybeans, wheat, corn, cotton, oats, beans). If urban glyphosate usage data is available, I plan to use graduated symbols to highlight cities where its use is highest.

TO-DO

– does “No data” mean 0 or that they didn’t look at those states? – Look at the AHS air conditioning data: https://www.census.gov/programs-surveys/ahs/data/interactive/ahstablecreator.html?s_areas=42660&s_year=2021&s_tablename=TABLE3&s_bygroup1=1&s_bygroup2=1&s_filtergroup1=1&s_filtergroup2=1 – Figure out how to outline each county – Get the actual best heat measure that I want to use – Most prevalent type of gly resistant crop for county – Urban gly usage data?

Functions

plot_gly <- function(year, est_type) {
  if (year==2013) {
    data = merged_2013
  }
  if (year==2014) {
    data = merged_2014
   }
  if (year==2015) {
    data = merged_2015
  }
  if (year==2016) {
    data = merged_2016
  }
  if (year==2017) {
    data = merged_2017
  }
  
  if (est_type=="low") {
    map <- ggplot(data) +
    geom_sf(aes(fill = gly_per_land_low), linetype=0) +
    ggtitle("Glyphosate per Land Area (Low Est.)", subtitle=year) +
    theme_void(base_size = 14) +
    theme(
      plot.title = element_text(hjust = 0.5),
      axis.text.x = element_blank(),
      axis.text.y = element_blank(),
      axis.ticks = element_blank(),
      rect = element_blank()
    ) +
    scale_fill_gradientn("Glyphosate usage \n(kg/km^2)",
                         colors = met.brewer("Tam"),
                         na.value = "grey50")
  }
  
  if (est_type=="high") {
    map <- ggplot(data) +
    geom_sf(aes(fill = gly_per_land_high), linetype=0) +
    ggtitle("Glyphosate per Land Area (High Est.)", subtitle=year) +
    theme_void(base_size = 14) +
    theme(
      plot.title = element_text(hjust = 0.5),
      axis.text.x = element_blank(),
      axis.text.y = element_blank(),
      axis.ticks = element_blank(),
      rect = element_blank()
    ) +
    scale_fill_gradientn("Glyphosate usage \n(kg/km^2)",
                         colors = met.brewer("Tam"),
                         na.value = "grey50")
  }
  map
}

compounddata_formap <- function(compound, year) {
  usgs_data <- usgs_ests %>% filter(COMPOUND==compound) %>% 
    filter(YEAR==year) %>% mutate(STATE_FIPS_CODE = 
                        str_pad(as.character(STATE_FIPS_CODE),
                              2, pad = "0"),
                COUNTY_FIPS_CODE =                  
                  str_pad(as.character(COUNTY_FIPS_CODE),
                              3, pad = "0"),
                fips = as.factor(paste0(STATE_FIPS_CODE, 
                                        COUNTY_FIPS_CODE)))
  
  if (year==2013) {
    us_county_data <- us_county_2013
    prism_county_data <- prism_county_2013_85s
    svi_data <- svi_2014
  }
  
  if (year==2014) {
    prism_county_data <- prism_county_2014_85s
    svi_data <- svi_2014
  }
  
  if (year==2015) {
    prism_county_data <- prism_county_2015_85s
    svi_data <- svi_2016
  }
  
  if (year==2016) {
    prism_county_data <- prism_county_2016_85s
    svi_data <- svi_2016
  }
  
  if (year==2017) {
    prism_county_data <- prism_county_2017_85s
    svi_data <- svi_2018
  }
  
  merged_data <- left_join(us_county_2013, prism_county_data) %>% 
  left_join(usgs_data) %>% filter(STATEFP != "02" & STATEFP != "15" &
                                   STATEFP != "72") %>%
  mutate(epest_per_land_low = EPEST_LOW_KG/(ALAND/1e+6),
         epest_per_land_high = EPEST_HIGH_KG/(ALAND/1e+6),
         epest_per_totalarea_low = EPEST_LOW_KG/((ALAND+AWATER)/1e+6),
         epest_per_totalarea_high = EPEST_HIGH_KG/((ALAND+AWATER)/1e+6)) %>%
  left_join(svi_data)

  merged_data
}

hot_compounddata_formap <- function(compound, year, temp_c) {
  compounddata_formap(compound, year) %>% filter(tmean85_spray>=temp_c)
}

all_leaflet_low <- function(compound, year) {

  #get data
  data <- compounddata_formap(compound, year)
  
  map <- data %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Voyager") %>%
    addPolygons(
             stroke = FALSE,
             smoothFactor = 0,
             fillOpacity = 0.5,
             fillColor= ~
               colorNumeric("YlOrRd", epest_per_land_low)(epest_per_land_low)) %>%
  addLegend(pal = colorNumeric("YlOrRd", data$epest_per_land_low),  # Define color palette
        values = data$epest_per_land_low,  # Values to map
        title = "Pesticide Usage (kg/km^2)",  # Legend title
        position = "bottomright",  # Position of the legend
        na.label = "No data"
    )
  
  map
}

all_leaflet_high <- function(compound, year) {

  #get data
  data <- compounddata_formap(compound, year)
  
  map <- data %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Voyager") %>%
    addPolygons(
             stroke = FALSE,
             smoothFactor = 0,
             fillOpacity = 0.5,
             fillColor= ~
               colorNumeric("YlOrRd", epest_per_land_high)(epest_per_land_high)) %>%
  addLegend(pal = colorNumeric("YlOrRd", data$epest_per_land_high),  # Define color palette
        values = data$epest_per_land_high,  # Values to map
        title = "Pesticide Usage (kg/km^2)",  # Legend title
        position = "bottomright",  # Position of the legend
        na.label = "No data"
    )
  
  map
}

hot_leaflet_low <- function(compound, year, temp_c) {

  #get data
  data <- hot_compounddata_formap(compound, year, temp_c)
  
  map <- data %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Voyager") %>%
    addPolygons(
             stroke = FALSE,
             smoothFactor = 0,
             fillOpacity = 0.5,
             fillColor= ~
               colorNumeric("YlOrRd", epest_per_land_low)(epest_per_land_low)) %>%
  addLegend(pal = colorNumeric("YlOrRd", data$epest_per_land_low),  # Define color palette
        values = data$epest_per_land_low,  # Values to map
        title = "Pesticide Usage (kg/km^2)",  # Legend title
        position = "bottomright",  # Position of the legend
        na.label = "No data"
    )
  
  map
}

hot_leaflet_high <- function(compound, year, temp_c) {

  #get data
  data <- hot_compounddata_formap(compound, year, temp_c)
  
  map <- data %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Voyager") %>%
    addPolygons(
             stroke = FALSE,
             smoothFactor = 0,
             fillOpacity = 0.5,
             fillColor= ~
               colorNumeric("YlOrRd", epest_per_land_high)(epest_per_land_high)) %>%
  addLegend(pal = colorNumeric("YlOrRd", data$epest_per_land_high),  # Define color palette
        values = data$epest_per_land_high,  # Values to map
        title = "Pesticide Usage (kg/km^2)",  # Legend title
        position = "bottomright",  # Position of the legend
        na.label = "No data"
    )
  
  map
}

Data Used

I am using USGS pesticide estimate data, Census Bureau data on U.S. county boundaries, and historical heat records in the U.S. from PRISM.

USGS Estimated Annual Agricultural Pesticide Use for Counties of the Conterminous United States, 2013-2017 (ver. 2.0, May 2020) data were obtained from here: https://www.sciencebase.gov/catalog/item/5e95c12282ce172707f2524e. Cartographic Boundary Files from the U.S Census Bureau were obtained from here for the year 2013: https://www.census.gov/geographies/mapping-files/time-series/geo/carto-boundary-file.2013.html#list-tab-1556094155. PRISM daily mean temperature data for 2013 was obtained from Robbie Parks’ Github (https://github.com/rmp15/PRISM-grids-into-FIPS-ZIP-censustract-USA/tree/main/output/fips/tmean). Citation for this work is Tuholske, C., Lynch, V.D., Spriggs, R. et al. Hazardous heat exposure among incarcerated people in the United States. Nat Sustain 7, 394–398 (2024). https://doi.org/10.1038/s41893-024-01293-y. Data dictionaries can be found at the above links.

I obtained CDC SVI data from: https://www.atsdr.cdc.gov/placeandhealth/svi/data_documentation_download.html. I obtained data for 2014, 2016, and 2018, and will use 2014 for years 2013-14, 2016 for 2015-16, and 2018 for 2017-2018. Data documentation is available at the same link (https://www.atsdr.cdc.gov/placeandhealth/svi/data_documentation_download.html).

#usgs_1992to2012 #read all these files from the list of files, then rbind them
usgs_ests <- fread(here("mapping data/EPest_county_estimates_2013_2017_v2.txt"))
usgs_compounds <- usgs_ests %>% select(c(COMPOUND)) %>% unique()
usgs_gly <- usgs_ests %>% filter(COMPOUND=="GLYPHOSATE")
usgs_gly_2013 <- usgs_gly %>% filter(YEAR==2013)
usgs_gly_2014 <- usgs_gly %>% filter(YEAR==2014)
usgs_gly_2015 <- usgs_gly %>% filter(YEAR==2015)
usgs_gly_2016 <- usgs_gly %>% filter(YEAR==2016)
usgs_gly_2017 <- usgs_gly %>% filter(YEAR==2017)

glimpse(usgs_gly)
## Rows: 15,314
## Columns: 6
## $ COMPOUND         <chr> "GLYPHOSATE", "GLYPHOSATE", "GLYPHOSATE", "GLYPHOSATE…
## $ YEAR             <int> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013,…
## $ STATE_FIPS_CODE  <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ COUNTY_FIPS_CODE <int> 1, 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29…
## $ EPEST_LOW_KG     <dbl> 13336.7, 48158.8, 9683.4, 379.4, 5589.1, 7321.7, 2719…
## $ EPEST_HIGH_KG    <dbl> 13538.5, 49014.8, 9683.4, 380.9, 5675.6, 7442.3, 2825…
dim(usgs_gly)
## [1] 15314     6
us_county_2013 <- st_read(here("mapping data/cb_2013_us_county_20m/cb_2013_us_county_20m.shp"))
## Reading layer `cb_2013_us_county_20m' from data source 
##   `/home/NETID/cm763/pesticide-heat-mapping/mapping data/cb_2013_us_county_20m/cb_2013_us_county_20m.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 3221 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -179.1473 ymin: 17.88481 xmax: 179.7785 ymax: 71.35256
## Geodetic CRS:  NAD83
glimpse(us_county_2013)
## Rows: 3,221
## Columns: 10
## $ STATEFP  <chr> "06", "13", "18", "21", "21", "26", "01", "48", "01", "36", "…
## $ COUNTYFP <chr> "049", "257", "127", "077", "053", "149", "085", "485", "037"…
## $ COUNTYNS <chr> "00277289", "00350028", "00450382", "00516885", "00516873", "…
## $ AFFGEOID <chr> "0500000US06049", "0500000US13257", "0500000US18127", "050000…
## $ GEOID    <chr> "06049", "13257", "18127", "21077", "21053", "26149", "01085"…
## $ NAME     <chr> "Modoc", "Stephens", "Porter", "Gallatin", "Clinton", "St. Jo…
## $ LSAD     <chr> "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "…
## $ ALAND    <dbl> 10140841147, 463948108, 1083011641, 254698298, 510864239, 129…
## $ AWATER   <dbl> 745929078, 13121937, 268393387, 16519758, 21164150, 52934828,…
## $ geometry <MULTIPOLYGON [°]> MULTIPOLYGON (((-121.4475 4..., MULTIPOLYGON (((…
dim(us_county_2013)
## [1] 3221   10
prism_county_2013 <- readRDS(here("mapping data/weighted_area_raster_fips_tmean_daily_2013.rds"))
glimpse(prism_county_2013)
## Rows: 1,134,420
## Columns: 6
## $ fips  <fct> 01001, 01003, 01005, 01007, 01009, 01011, 01013, 01015, 01017, 0…
## $ tmean <dbl> 6.188720, 9.711895, 7.783865, 6.224829, 5.372772, 6.969209, 7.23…
## $ date  <chr> "01/01/2013", "01/01/2013", "01/01/2013", "01/01/2013", "01/01/2…
## $ day   <chr> "01", "01", "01", "01", "01", "01", "01", "01", "01", "01", "01"…
## $ month <chr> "01", "01", "01", "01", "01", "01", "01", "01", "01", "01", "01"…
## $ year  <dbl> 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013…
dim(prism_county_2013)
## [1] 1134420       6
prism_county_2014 <- readRDS(here("mapping data/weighted_area_raster_fips_tmean_daily_2014.rds"))

prism_county_2015 <- readRDS(here("mapping data/weighted_area_raster_fips_tmean_daily_2015.rds"))

prism_county_2016 <- readRDS(here("mapping data/weighted_area_raster_fips_tmean_daily_2016.rds"))

prism_county_2017 <- readRDS(here("mapping data/weighted_area_raster_fips_tmean_daily_2017.rds"))

svi_2014 <- fread(here("mapping data/SVI_2014_US_county.csv"))
glimpse(svi_2014)
## Rows: 3,142
## Columns: 127
## $ AFFGEOID           <chr> "0500000US01001", "0500000US01003", "0500000US01005…
## $ ST                 <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ STATE              <chr> "ALABAMA", "ALABAMA", "ALABAMA", "ALABAMA", "ALABAM…
## $ ST_ABBR            <chr> "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL…
## $ COUNTY             <chr> "Autauga", "Baldwin", "Barbour", "Bibb", "Blount", …
## $ FIPS               <int> 1001, 1003, 1005, 1007, 1009, 1011, 1013, 1015, 101…
## $ LOCATION           <chr> "Autauga County, Alabama", "Baldwin County, Alabama…
## $ AREA_SQMI          <dbl> 594.4366, 1589.8073, 884.8767, 622.5824, 644.8065, …
## $ E_TOTPOP           <dbl> 55136, 191205, 27119, 22653, 57645, 10693, 20523, 1…
## $ M_TOTPOP           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ E_HU               <dbl> 22431, 105563, 11833, 8985, 23868, 4469, 9934, 5330…
## $ M_HU               <dbl> 67, 168, 120, 66, 77, 100, 67, 209, 50, 89, 68, 33,…
## $ E_HH               <dbl> 20304, 73058, 9145, 7078, 20934, 3746, 8253, 45348,…
## $ M_HH               <dbl> 458, 1241, 311, 390, 399, 216, 252, 705, 393, 451, …
## $ E_POV              <dbl> 7006, 25988, 5832, 3596, 9866, 2085, 5239, 24794, 8…
## $ M_POV              <dbl> 935, 2457, 639, 770, 947, 477, 517, 1491, 823, 876,…
## $ E_UNEMP            <dbl> 2252, 7856, 1527, 975, 2291, 809, 1075, 7257, 1916,…
## $ M_UNEMP            <dbl> 352, 918, 282, 310, 358, 187, 177, 594, 304, 221, 3…
## $ E_PCI              <dbl> 24644, 26851, 17350, 18110, 20501, 17706, 18115, 21…
## $ M_PCI              <dbl> 780, 712, 821, 1477, 719, 1557, 832, 573, 1482, 138…
## $ E_NOHSDP           <dbl> 5012, 14615, 4790, 3466, 8567, 2511, 3288, 15674, 5…
## $ M_NOHSDP           <dbl> 561, 1051, 341, 500, 697, 319, 278, 787, 388, 373, …
## $ E_AGE65            <dbl> 7321, 33782, 4180, 3209, 9172, 1529, 3592, 17820, 6…
## $ M_AGE65            <dbl> 84, 151, 81, 84, 89, 0, 72, 95, 38, 94, 135, 0, 0, …
## $ E_AGE17            <dbl> 14214, 43186, 5862, 4903, 13877, 2283, 4835, 26337,…
## $ M_AGE17            <dbl> 0, 0, 27, 0, 33, 0, 50, 0, 26, 78, 50, 0, 0, 46, 10…
## $ E_DISABL           <dbl> 8700, 26603, 4856, 3343, 9729, 2049, 4301, 21597, 6…
## $ M_DISABL           <dbl> 674, 1468, 379, 429, 748, 309, 284, 861, 507, 552, …
## $ E_SNGPNT           <dbl> 1649, 5027, 1098, 637, 1650, 535, 1104, 4406, 1725,…
## $ M_SNGPNT           <dbl> 300.2, 603.7, 160.7, 211.3, 299.2, 179.6, 173.7, 50…
## $ E_MINRTY           <dbl> 13103, 32158, 14614, 5717, 6825, 8321, 9517, 31315,…
## $ M_MINRTY           <dbl> 109, 162, 166, 21, 52, 18, 22, 68, 24, 21, 91, 18, …
## $ E_LIMENG           <dbl> 249, 2571, 549, 112, 954, 257, 82, 915, 67, 25, 113…
## $ M_LIMENG           <dbl> 138.5, 658.6, 261.6, 110.9, 222.2, 250.7, 100.2, 24…
## $ E_MUNIT            <dbl> 842, 18988, 129, 163, 215, 65, 136, 2025, 755, 195,…
## $ M_MUNIT            <dbl> 344.4, 1015.0, 68.6, 97.6, 95.2, 50.6, 67.7, 296.9,…
## $ E_MOBILE           <dbl> 4401, 12200, 3175, 2379, 5497, 1438, 2397, 7467, 30…
## $ M_MOBILE           <dbl> 367, 821, 270, 314, 482, 228, 203, 511, 320, 398, 4…
## $ E_CROWD            <dbl> 530, 998, 178, 17, 344, 122, 95, 575, 506, 266, 567…
## $ M_CROWD            <dbl> 231.2, 265.6, 83.6, 25.8, 142.4, 118.7, 49.2, 165.0…
## $ E_NOVEH            <dbl> 1081, 2242, 802, 299, 823, 703, 701, 2562, 1238, 34…
## $ M_NOVEH            <dbl> 218, 308, 152, 103, 180, 188, 127, 307, 220, 98, 18…
## $ E_GROUPQ           <dbl> 442, 2611, 2869, 1576, 572, 576, 270, 2768, 447, 39…
## $ M_GROUPQ           <dbl> 133, 405, 290, 177, 147, 248, 146, 500, 105, 187, 1…
## $ EP_POV             <dbl> 12.8, 13.8, 24.1, 17.0, 17.3, 20.5, 25.9, 21.7, 23.…
## $ MP_POV             <dbl> 1.7, 1.3, 2.6, 3.6, 1.7, 4.7, 2.6, 1.3, 2.4, 3.4, 2…
## $ EP_UNEMP           <dbl> 8.5, 8.6, 14.2, 10.9, 9.3, 17.4, 12.2, 13.5, 12.7, …
## $ MP_UNEMP           <dbl> 1.3, 1.0, 2.5, 3.4, 1.4, 3.6, 1.9, 1.1, 2.0, 2.0, 1…
## $ EP_PCI             <dbl> 24644, 26851, 17350, 18110, 20501, 17706, 18115, 21…
## $ MP_PCI             <dbl> 780, 712, 821, 1477, 719, 1557, 832, 573, 1482, 138…
## $ EP_NOHSDP          <dbl> 13.8, 11.0, 25.4, 22.1, 21.9, 34.5, 23.4, 19.9, 22.…
## $ MP_NOHSDP          <dbl> 1.5, 0.8, 1.8, 3.2, 1.8, 4.4, 2.0, 1.0, 1.6, 2.0, 1…
## $ EP_AGE65           <dbl> 13.3, 17.7, 15.4, 14.2, 15.9, 14.3, 17.5, 15.2, 17.…
## $ MP_AGE65           <dbl> 0.2, 0.1, 0.3, 0.4, 0.2, 0.0, 0.3, 0.1, 0.1, 0.4, 0…
## $ EP_AGE17           <dbl> 25.8, 22.6, 21.6, 21.6, 24.1, 21.4, 23.6, 22.5, 21.…
## $ MP_AGE17           <dbl> 0.0, 0.0, 0.1, 0.0, 0.1, 0.0, 0.2, 0.0, 0.1, 0.3, 0…
## $ EP_DISABL          <dbl> 16.0, 14.1, 20.0, 15.8, 17.0, 20.2, 21.2, 18.7, 19.…
## $ MP_DISABL          <dbl> 1.2, 0.8, 1.5, 2.0, 1.3, 3.0, 1.4, 0.7, 1.5, 2.2, 1…
## $ EP_SNGPNT          <dbl> 8.1, 6.9, 12.0, 9.0, 7.9, 14.3, 13.4, 9.7, 12.4, 8.…
## $ MP_SNGPNT          <dbl> 1.5, 0.8, 1.7, 2.9, 1.4, 4.7, 2.1, 1.1, 1.6, 2.0, 1…
## $ EP_MINRTY          <dbl> 23.8, 16.8, 53.9, 25.2, 11.8, 77.8, 46.4, 26.7, 42.…
## $ MP_MINRTY          <dbl> 0.2, 0.1, 0.6, 0.1, 0.1, 0.2, 0.1, 0.1, 0.1, 0.1, 0…
## $ EP_LIMENG          <dbl> 0.5, 1.4, 2.2, 0.5, 1.8, 2.6, 0.4, 0.8, 0.2, 0.1, 2…
## $ MP_LIMENG          <dbl> 0.3, 0.4, 1.0, 0.5, 0.4, 2.5, 0.5, 0.2, 0.3, 0.3, 0…
## $ EP_MUNIT           <dbl> 3.8, 18.0, 1.1, 1.8, 0.9, 1.5, 1.4, 3.8, 4.5, 1.2, …
## $ MP_MUNIT           <dbl> 1.5, 1.0, 0.6, 1.1, 0.4, 1.1, 0.7, 0.6, 1.0, 0.5, 0…
## $ EP_MOBILE          <dbl> 19.6, 11.6, 26.8, 26.5, 23.0, 32.2, 24.1, 14.0, 18.…
## $ MP_MOBILE          <dbl> 1.6, 0.8, 2.3, 3.5, 2.0, 5.0, 2.1, 1.0, 1.9, 2.5, 2…
## $ EP_CROWD           <dbl> 2.6, 1.4, 1.9, 0.2, 1.6, 3.3, 1.2, 1.3, 3.6, 2.3, 3…
## $ MP_CROWD           <dbl> 1.1, 0.4, 0.9, 0.4, 0.7, 3.2, 0.6, 0.4, 1.0, 1.2, 1…
## $ EP_NOVEH           <dbl> 5.3, 3.1, 8.8, 4.2, 3.9, 18.8, 8.5, 5.6, 8.9, 3.0, …
## $ MP_NOVEH           <dbl> 1.1, 0.4, 1.6, 1.5, 0.9, 5.0, 1.5, 0.7, 1.6, 0.8, 1…
## $ EP_GROUPQ          <dbl> 0.8, 1.4, 10.6, 7.0, 1.0, 5.4, 1.3, 2.4, 1.3, 1.5, …
## $ MP_GROUPQ          <dbl> 0.2, 0.2, 1.1, 0.8, 0.3, 2.3, 0.7, 0.4, 0.3, 0.7, 0…
## $ EPL_POV            <dbl> 0.2833, 0.3531, 0.8736, 0.5565, 0.5782, 0.7561, 0.9…
## $ EPL_UNEMP          <dbl> 0.5174, 0.5288, 0.9319, 0.7650, 0.6074, 0.9787, 0.8…
## $ EPL_PCI            <dbl> 0.3913, 0.2467, 0.9274, 0.8949, 0.7348, 0.9115, 0.8…
## $ EPL_NOHSDP         <dbl> 0.5126, 0.3225, 0.9236, 0.8360, 0.8287, 0.9901, 0.8…
## $ SPL_THEME1         <dbl> 1.7046, 1.4511, 3.6565, 3.0525, 2.7491, 3.6364, 3.5…
## $ RPL_THEME1         <dbl> 0.4145, 0.3372, 0.9628, 0.8201, 0.7348, 0.9580, 0.9…
## $ EPL_AGE65          <dbl> 0.1980, 0.6281, 0.3964, 0.2677, 0.4429, 0.2767, 0.6…
## $ EPL_AGE17          <dbl> 0.8532, 0.4597, 0.3327, 0.3368, 0.6714, 0.3044, 0.6…
## $ EPL_DISABL         <dbl> 0.5762, 0.3983, 0.8475, 0.5571, 0.6587, 0.8580, 0.8…
## $ EPL_SNGPNT         <dbl> 0.4330, 0.2372, 0.8937, 0.5804, 0.3903, 0.9659, 0.9…
## $ SPL_THEME2         <dbl> 2.0605, 1.7233, 2.4702, 1.7421, 2.1633, 2.4050, 3.0…
## $ RPL_THEME2         <dbl> 0.5387, 0.2811, 0.8195, 0.2948, 0.6199, 0.7826, 0.9…
## $ EPL_MINRTY         <dbl> 0.6332, 0.5336, 0.9086, 0.6517, 0.4247, 0.9819, 0.8…
## $ EPL_LIMENG         <dbl> 0.3722, 0.6778, 0.7720, 0.3961, 0.7297, 0.8042, 0.3…
## $ SPL_THEME3         <dbl> 1.0054, 1.2114, 1.6807, 1.0478, 1.1544, 1.7861, 1.2…
## $ RPL_THEME3         <dbl> 0.4986, 0.6288, 0.8898, 0.5237, 0.5864, 0.9322, 0.6…
## $ EPL_MUNIT          <dbl> 0.6135, 0.9733, 0.2006, 0.3499, 0.1649, 0.2767, 0.2…
## $ EPL_MOBILE         <dbl> 0.7701, 0.5183, 0.9064, 0.9019, 0.8488, 0.9577, 0.8…
## $ EPL_CROWD          <dbl> 0.7011, 0.3098, 0.5317, 0.0216, 0.4167, 0.8103, 0.2…
## $ EPL_NOVEH          <dbl> 0.3820, 0.0844, 0.8341, 0.2073, 0.1713, 0.9898, 0.8…
## $ EPL_GROUPQ         <dbl> 0.0863, 0.2865, 0.9414, 0.8781, 0.1407, 0.8259, 0.2…
## $ SPL_THEME4         <dbl> 2.5530, 2.1722, 3.4142, 2.3588, 1.7424, 3.8602, 2.4…
## $ RPL_THEME4         <dbl> 0.5100, 0.3238, 0.9089, 0.4152, 0.1563, 0.9866, 0.4…
## $ SPL_THEMES         <dbl> 7.3235, 6.5581, 11.2216, 8.2012, 7.8093, 11.6877, 1…
## $ RPL_THEMES         <dbl> 0.4696, 0.3432, 0.9742, 0.6278, 0.5581, 0.9914, 0.9…
## $ F_POV              <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, …
## $ F_UNEMP            <dbl> 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, …
## $ F_PCI              <dbl> 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ F_NOHSDP           <dbl> 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, …
## $ F_THEME1           <dbl> 0, 0, 3, 0, 0, 3, 1, 1, 0, 0, 0, 2, 2, 1, 1, 0, 0, …
## $ F_AGE65            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ F_AGE17            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ F_DISABL           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, …
## $ F_SNGPNT           <dbl> 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ F_THEME2           <dbl> 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, …
## $ F_MINRTY           <dbl> 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ F_LIMENG           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ F_THEME3           <dbl> 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ F_MUNIT            <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ F_MOBILE           <dbl> 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, 1, 0, 0, …
## $ F_CROWD            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ F_NOVEH            <dbl> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, …
## $ F_GROUPQ           <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ F_THEME4           <dbl> 0, 1, 2, 1, 0, 2, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, …
## $ F_TOTAL            <dbl> 0, 1, 6, 1, 0, 7, 2, 1, 1, 2, 1, 4, 3, 2, 2, 0, 0, …
## $ E_UNINSUR          <dbl> 5814, 27758, 3760, 2725, 6922, 1791, 2603, 15258, 5…
## $ M_UNINSUR          <dbl> 674, 2097, 412, 554, 660, 313, 327, 1031, 476, 504,…
## $ EP_UNINSUR         <dbl> 10.7, 14.7, 15.5, 12.9, 12.1, 17.6, 12.8, 13.2, 15.…
## $ MP_UNINSUR         <dbl> 1.2, 1.1, 1.7, 2.6, 1.2, 3.0, 1.6, 0.9, 1.4, 2.0, 1…
## $ E_DAYPOP           <dbl> 43534, 177010, 29769, 19274, 41857, 10639, 19896, 1…
## $ Shape              <chr> "(-86.64274304015402, 32.53492165233667)", "(-87.72…
## $ `Shape.STArea()`   <dbl> 0.1502601, 0.4099377, 0.2232557, 0.1565251, 0.16440…
## $ `Shape.STLength()` <dbl> 2.052620, 4.277985, 2.566611, 1.886957, 2.392305, 2…
dim(svi_2014)
## [1] 3142  127
svi_2016 <- fread(here("mapping data/SVI_2016_US_county.csv"))
glimpse(svi_2016)
## Rows: 3,142
## Columns: 123
## $ ST         <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ STATE      <chr> "ALABAMA", "ALABAMA", "ALABAMA", "ALABAMA", "ALABAMA", "ALA…
## $ ST_ABBR    <chr> "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL",…
## $ COUNTY     <chr> "Autauga", "Blount", "Chambers", "Coffee", "Colbert", "Covi…
## $ FIPS       <int> 1001, 1009, 1017, 1031, 1033, 1039, 1043, 1045, 1051, 1055,…
## $ LOCATION   <chr> "Autauga County, Alabama", "Blount County, Alabama", "Chamb…
## $ AREA_SQMI  <dbl> 594.4461, 644.8065, 596.5311, 678.9857, 592.6197, 1030.4558…
## $ E_TOTPOP   <dbl> 55049, 57704, 34018, 50991, 54377, 37729, 81316, 49607, 809…
## $ M_TOTPOP   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ E_HU       <dbl> 22714, 23850, 16905, 22862, 26156, 18783, 37150, 22839, 333…
## $ M_HU       <dbl> 75, 59, 53, 106, 88, 66, 92, 95, 217, 126, 63, 138, 63, 146…
## $ E_HH       <dbl> 20800, 20619, 13851, 19375, 22105, 15179, 31081, 18794, 288…
## $ M_HH       <dbl> 391, 403, 359, 446, 428, 337, 517, 416, 541, 505, 246, 369,…
## $ E_POV      <dbl> 6697, 9441, 6805, 8219, 8910, 7031, 13718, 9897, 10893, 187…
## $ M_POV      <dbl> 1037, 963, 849, 901, 839, 644, 1340, 914, 1217, 1295, 455, …
## $ E_UNEMP    <dbl> 1437, 1367, 1136, 1410, 1795, 1723, 2081, 2145, 2694, 3837,…
## $ M_UNEMP    <dbl> 277, 284, 232, 235, 301, 314, 314, 302, 522, 462, 170, 224,…
## $ E_PCI      <dbl> 26168, 21033, 21532, 25325, 23318, 21738, 21041, 22834, 247…
## $ M_PCI      <dbl> 1221, 689, 1399, 939, 854, 826, 734, 783, 1013, 489, 1459, …
## $ E_NOHSDP   <dbl> 4528, 7882, 4708, 5145, 6344, 5051, 10050, 4726, 7238, 1258…
## $ M_NOHSDP   <dbl> 445, 645, 421, 418, 506, 345, 745, 396, 594, 619, 239, 408,…
## $ E_AGE65    <dbl> 7695, 9921, 6255, 8048, 10034, 7613, 14435, 7654, 11243, 17…
## $ M_AGE65    <dbl> 104, 123, 48, 89, 71, 69, 106, 79, 82, 75, 46, 0, 29, 78, 5…
## $ E_AGE17    <dbl> 13853, 13601, 7283, 12122, 11735, 8318, 18279, 11728, 18553…
## $ M_AGE17    <dbl> 34, 29, 52, 0, 0, 67, 110, 0, 0, 107, 47, 0, 66, 24, 60, 63…
## $ E_DISABL   <dbl> 10009, 8538, 5960, 8942, 10561, 7858, 12621, 9865, 13201, 2…
## $ M_DISABL   <dbl> 850, 663, 492, 597, 677, 439, 739, 494, 967, 937, 296, 462,…
## $ E_SNGPNT   <dbl> 1516, 1614, 1354, 2018, 1879, 1240, 2168, 1649, 2688, 3631,…
## $ M_SNGPNT   <dbl> 267.0, 301.0, 187.6, 266.6, 290.0, 171.4, 337.1, 251.6, 394…
## $ E_MINRTY   <dbl> 13386, 7122, 14715, 14598, 11499, 6220, 6128, 14787, 21301,…
## $ M_MINRTY   <dbl> 161, 147, 24, 7, 46, 6, 104, 77, 109, 108, 18, 55, 91, 58, …
## $ E_LIMENG   <dbl> 432, 1018, 114, 716, 104, 211, 746, 395, 635, 1475, 115, 73…
## $ M_LIMENG   <dbl> 163.3, 248.1, 100.7, 242.1, 122.2, 112.3, 207.7, 183.9, 247…
## $ E_MUNIT    <dbl> 1034, 190, 597, 218, 641, 174, 901, 304, 894, 1515, 86, 199…
## $ M_MUNIT    <dbl> 329.9, 85.8, 170.3, 81.8, 154.7, 75.9, 191.1, 113.8, 235.8,…
## $ E_MOBILE   <dbl> 4095, 5467, 2695, 2863, 2478, 3505, 7996, 4020, 5587, 5487,…
## $ M_MOBILE   <dbl> 379, 429, 317, 202, 333, 255, 575, 302, 472, 441, 239, 297,…
## $ E_CROWD    <dbl> 254, 391, 555, 270, 150, 237, 547, 269, 332, 639, 59, 479, …
## $ M_CROWD    <dbl> 104.5, 130.6, 150.8, 101.0, 76.9, 90.7, 136.0, 101.2, 129.8…
## $ E_NOVEH    <dbl> 1024, 816, 1110, 1126, 1361, 797, 1264, 1198, 1291, 2513, 4…
## $ M_NOVEH    <dbl> 242, 198, 183, 204, 224, 123, 218, 213, 308, 286, 117, 155,…
## $ E_GROUPQ   <dbl> 490, 552, 482, 616, 432, 612, 1164, 1209, 4651, 1403, 268, …
## $ M_GROUPQ   <dbl> 163, 131, 106, 134, 102, 151, 216, 219, 570, 249, 79, 170, …
## $ EP_POV     <dbl> 12.3, 16.5, 20.3, 16.4, 16.5, 19.0, 17.1, 20.5, 14.3, 18.4,…
## $ MP_POV     <dbl> 1.9, 1.7, 2.5, 1.8, 1.6, 1.7, 1.7, 1.9, 1.6, 1.3, 2.7, 0.7,…
## $ EP_UNEMP   <dbl> 5.6, 6.0, 7.5, 6.2, 7.5, 10.7, 6.1, 10.4, 7.5, 8.5, 7.5, 8.…
## $ MP_UNEMP   <dbl> 1.1, 1.2, 1.5, 1.0, 1.3, 1.9, 0.9, 1.4, 1.4, 1.0, 2.2, 0.5,…
## $ EP_PCI     <dbl> 26168, 21033, 21532, 25325, 23318, 21738, 21041, 22834, 247…
## $ MP_PCI     <dbl> 1221, 689, 1399, 939, 854, 826, 734, 783, 1013, 489, 1459, …
## $ EP_NOHSDP  <dbl> 12.4, 20.0, 19.7, 14.9, 16.6, 19.1, 17.8, 14.2, 13.2, 17.5,…
## $ MP_NOHSDP  <dbl> 1.2, 1.6, 1.8, 1.2, 1.3, 1.3, 1.3, 1.2, 1.1, 0.9, 1.9, 0.6,…
## $ EP_AGE65   <dbl> 14.0, 17.2, 18.4, 15.8, 18.5, 20.2, 17.8, 15.4, 13.9, 17.4,…
## $ MP_AGE65   <dbl> 0.2, 0.2, 0.1, 0.2, 0.1, 0.2, 0.1, 0.2, 0.1, 0.1, 0.3, 0.0,…
## $ EP_AGE17   <dbl> 25.2, 23.6, 21.4, 23.8, 21.6, 22.0, 22.5, 23.6, 22.9, 22.0,…
## $ MP_AGE17   <dbl> 0.1, 0.1, 0.2, 0.0, 0.0, 0.2, 0.1, 0.0, 0.0, 0.1, 0.3, 0.0,…
## $ EP_DISABL  <dbl> 18.4, 14.9, 17.7, 18.2, 19.6, 21.1, 15.7, 21.0, 17.4, 19.8,…
## $ MP_DISABL  <dbl> 1.6, 1.2, 1.5, 1.2, 1.3, 1.2, 0.9, 1.1, 1.3, 0.9, 1.7, 0.4,…
## $ EP_SNGPNT  <dbl> 7.3, 7.8, 9.8, 10.4, 8.5, 8.2, 7.0, 8.8, 9.3, 9.2, 7.2, 11.…
## $ MP_SNGPNT  <dbl> 1.3, 1.5, 1.3, 1.4, 1.3, 1.1, 1.1, 1.3, 1.4, 0.9, 1.9, 0.6,…
## $ EP_MINRTY  <dbl> 24.3, 12.3, 43.3, 28.6, 21.1, 16.5, 7.5, 29.8, 26.3, 21.6, …
## $ MP_MINRTY  <dbl> 0.3, 0.3, 0.1, 0.0, 0.1, 0.0, 0.1, 0.2, 0.1, 0.1, 0.1, 0.1,…
## $ EP_LIMENG  <dbl> 0.8, 1.9, 0.4, 1.5, 0.2, 0.6, 1.0, 0.9, 0.8, 1.5, 0.7, 0.8,…
## $ MP_LIMENG  <dbl> 0.3, 0.5, 0.3, 0.5, 0.2, 0.3, 0.3, 0.4, 0.3, 0.3, 0.5, 0.1,…
## $ EP_MUNIT   <dbl> 4.6, 0.8, 3.5, 1.0, 2.5, 0.9, 2.4, 1.3, 2.7, 3.2, 1.0, 4.3,…
## $ MP_MUNIT   <dbl> 1.5, 0.4, 1.0, 0.4, 0.6, 0.4, 0.5, 0.5, 0.7, 0.4, 0.5, 0.4,…
## $ EP_MOBILE  <dbl> 18.0, 22.9, 15.9, 12.5, 9.5, 18.7, 21.5, 17.6, 16.8, 11.6, …
## $ MP_MOBILE  <dbl> 1.7, 1.8, 1.9, 0.9, 1.3, 1.4, 1.5, 1.3, 1.4, 0.9, 2.7, 0.6,…
## $ EP_CROWD   <dbl> 1.2, 1.9, 4.0, 1.4, 0.7, 1.6, 1.8, 1.4, 1.1, 1.6, 0.9, 1.2,…
## $ MP_CROWD   <dbl> 0.5, 0.6, 1.1, 0.5, 0.3, 0.6, 0.4, 0.5, 0.4, 0.5, 0.5, 0.2,…
## $ EP_NOVEH   <dbl> 4.9, 4.0, 8.0, 5.8, 6.2, 5.3, 4.1, 6.4, 4.5, 6.4, 6.7, 6.3,…
## $ MP_NOVEH   <dbl> 1.1, 1.0, 1.3, 1.0, 1.0, 0.8, 0.7, 1.1, 1.1, 0.7, 1.7, 0.4,…
## $ EP_GROUPQ  <dbl> 0.9, 1.0, 1.4, 1.2, 0.8, 1.6, 1.4, 2.4, 5.7, 1.4, 1.6, 0.9,…
## $ MP_GROUPQ  <dbl> 0.3, 0.2, 0.3, 0.3, 0.2, 0.4, 0.3, 0.4, 0.7, 0.2, 0.5, 0.2,…
## $ EPL_POV    <dbl> 0.2824, 0.5536, 0.7644, 0.5466, 0.5536, 0.7042, 0.5925, 0.7…
## $ EPL_UNEMP  <dbl> 0.3298, 0.3792, 0.5947, 0.4094, 0.5947, 0.8854, 0.3964, 0.8…
## $ EPL_PCI    <dbl> 0.3607, 0.7504, 0.7122, 0.4187, 0.5766, 0.6969, 0.7494, 0.6…
## $ EPL_NOHSDP <dbl> 0.4744, 0.8058, 0.7934, 0.6055, 0.6829, 0.7720, 0.7329, 0.5…
## $ SPL_THEME1 <dbl> 1.4473, 2.4890, 2.8647, 1.9803, 2.4078, 3.0586, 2.4712, 2.8…
## $ RPL_THEME1 <dbl> 0.3298, 0.6609, 0.7689, 0.5072, 0.6345, 0.8214, 0.6558, 0.7…
## $ EPL_AGE65  <dbl> 0.1964, 0.4909, 0.6151, 0.3464, 0.6243, 0.7676, 0.5524, 0.3…
## $ EPL_AGE17  <dbl> 0.8313, 0.6466, 0.3467, 0.6756, 0.3664, 0.4263, 0.4938, 0.6…
## $ EPL_DISABL <dbl> 0.7380, 0.4527, 0.6890, 0.7233, 0.8032, 0.8768, 0.5288, 0.8…
## $ EPL_SNGPNT <dbl> 0.3200, 0.4018, 0.7138, 0.7845, 0.5244, 0.4629, 0.2744, 0.5…
## $ SPL_THEME2 <dbl> 2.0856, 1.9920, 2.3645, 2.5298, 2.3184, 2.5336, 1.8494, 2.4…
## $ RPL_THEME2 <dbl> 0.5568, 0.4664, 0.7615, 0.8507, 0.7326, 0.8516, 0.3613, 0.7…
## $ EPL_MINRTY <dbl> 0.6339, 0.4238, 0.8402, 0.6877, 0.5880, 0.5205, 0.2786, 0.6…
## $ EPL_LIMENG <dbl> 0.5355, 0.7482, 0.2996, 0.6931, 0.1821, 0.4323, 0.5813, 0.5…
## $ SPL_THEME3 <dbl> 1.1694, 1.1719, 1.1398, 1.3808, 0.7701, 0.9529, 0.8599, 1.2…
## $ RPL_THEME3 <dbl> 0.5976, 0.5992, 0.5798, 0.7313, 0.3642, 0.4686, 0.4174, 0.6…
## $ EPL_MUNIT  <dbl> 0.6791, 0.1344, 0.5845, 0.1668, 0.4511, 0.1576, 0.4495, 0.2…
## $ EPL_MOBILE <dbl> 0.7268, 0.8465, 0.6734, 0.5613, 0.4368, 0.7488, 0.8138, 0.7…
## $ EPL_CROWD  <dbl> 0.2477, 0.5056, 0.8812, 0.3091, 0.0780, 0.3798, 0.4550, 0.3…
## $ EPL_NOVEH  <dbl> 0.3298, 0.1907, 0.7740, 0.4877, 0.5495, 0.3986, 0.2038, 0.5…
## $ EPL_GROUPQ <dbl> 0.1251, 0.1515, 0.3251, 0.2340, 0.0981, 0.4024, 0.3308, 0.5…
## $ SPL_THEME4 <dbl> 2.1086, 1.8287, 3.2381, 1.7590, 1.6135, 2.0872, 2.2528, 2.4…
## $ RPL_THEME4 <dbl> 0.2881, 0.1827, 0.8596, 0.1624, 0.1175, 0.2802, 0.3531, 0.4…
## $ SPL_THEMES <dbl> 6.8109, 7.4817, 9.6071, 7.6498, 7.1098, 8.6323, 7.4333, 8.9…
## $ RPL_THEMES <dbl> 0.3773, 0.4986, 0.8408, 0.5266, 0.4333, 0.6950, 0.4909, 0.7…
## $ F_POV      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_UNEMP    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_PCI      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_NOHSDP   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_THEME1   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_AGE65    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_AGE17    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_DISABL   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_SNGPNT   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_THEME2   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_MINRTY   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_LIMENG   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_THEME3   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_MUNIT    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_MOBILE   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_CROWD    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_NOVEH    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_GROUPQ   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_THEME4   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_TOTAL    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ E_UNINSUR  <dbl> 4852, 6388, 3979, 5253, 4932, 4617, 11704, 6130, 6665, 1236…
## $ M_UNINSUR  <dbl> 649, 740, 544, 464, 458, 448, 972, 562, 803, 857, 271, 597,…
## $ EP_UNINSUR <dbl> 8.9, 11.2, 11.8, 10.7, 9.1, 12.4, 14.5, 13.1, 8.8, 12.1, 11…
## $ MP_UNINSUR <dbl> 1.2, 1.3, 1.6, 0.9, 0.8, 1.2, 1.2, 1.2, 1.0, 0.8, 1.6, 0.6,…
## $ E_DAYPOP   <dbl> 40854, 42597, 27940, 47236, 56227, 37246, 79900, 49165, 669…
dim(svi_2016)
## [1] 3142  123
svi_2018 <- fread(here("mapping data/SVI_2018_US_county.csv"))
glimpse(svi_2018)
## Rows: 3,142
## Columns: 123
## $ ST         <int> 35, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ STATE      <chr> "NEW MEXICO", "ALABAMA", "ALABAMA", "ALABAMA", "ALABAMA", "…
## $ ST_ABBR    <chr> "NM", "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL",…
## $ COUNTY     <chr> "Rio Arriba", "Autauga", "Blount", "Butler", "Calhoun", "Ch…
## $ FIPS       <int> 35039, 1001, 1009, 1013, 1015, 1017, 1031, 1033, 1039, 1043…
## $ LOCATION   <chr> "Rio Arriba County, New Mexico", "Autauga County, Alabama",…
## $ AREA_SQMI  <dbl> 5860.8692, 594.4435, 644.8305, 776.8382, 605.8673, 596.5606…
## $ E_TOTPOP   <int> 39307, 55200, 57645, 20025, 115098, 33826, 51288, 54495, 37…
## $ M_TOTPOP   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ E_HU       <int> 20044, 23315, 24222, 10026, 53682, 16981, 23088, 26466, 189…
## $ M_HU       <int> 71, 71, 55, 51, 184, 72, 57, 58, 85, 92, 110, 192, 75, 44, …
## $ E_HH       <int> 12398, 21115, 20600, 6708, 45033, 13516, 19789, 21799, 1500…
## $ M_HH       <int> 439, 383, 396, 274, 683, 372, 339, 445, 381, 592, 402, 638,…
## $ E_POV      <int> -999, 8422, 8220, 4640, 20819, 5531, 7604, 8785, 6702, 1283…
## $ M_POV      <int> -999, 1137, 992, 521, 1317, 682, 731, 848, 683, 1294, 859, …
## $ E_UNEMP    <int> -999, 1065, 909, 567, 4628, 773, 1371, 1397, 1428, 1903, 17…
## $ M_UNEMP    <int> -999, 257, 193, 147, 526, 173, 244, 286, 252, 353, 275, 442…
## $ E_PCI      <int> -999, 29372, 22656, 20430, 24706, 22827, 27577, 24918, 2307…
## $ M_PCI      <int> -999, 2306, 905, 1258, 758, 1482, 976, 1055, 854, 906, 925,…
## $ E_NOHSDP   <int> 3669, 4204, 7861, 2141, 12620, 4383, 4833, 6016, 4530, 1040…
## $ M_NOHSDP   <int> 426, 475, 727, 268, 766, 356, 429, 555, 426, 776, 464, 675,…
## $ E_AGE65    <int> 7083, 8050, 10233, 3806, 19386, 6409, 8375, 10502, 7650, 14…
## $ M_AGE65    <int> 25, 75, 91, 21, 119, 85, 75, 73, 100, 95, 85, 111, 99, 36, …
## $ E_AGE17    <int> 9318, 13369, 13468, 4566, 25196, 7006, 12161, 11608, 8190, …
## $ M_AGE17    <int> 14, 32, 53, 88, 66, 94, 0, 38, 85, 140, 39, 0, 92, 0, 0, 18…
## $ E_DISABL   <int> 6280, 10465, 8114, 3492, 23598, 5570, 8788, 10086, 7703, 14…
## $ M_DISABL   <int> 495, 729, 592, 370, 1086, 422, 488, 520, 416, 958, 664, 860…
## $ E_SNGPNT   <int> 1330, 1586, 1437, 704, 4701, 1307, 1928, 1791, 1358, 2057, …
## $ M_SNGPNT   <dbl> 285.0, 319.9, 267.2, 143.9, 464.0, 244.6, 272.8, 278.2, 195…
## $ E_MINRTY   <int> 34397, 13788, 7413, 9641, 31675, 14954, 15065, 11630, 6205,…
## $ M_MINRTY   <dbl> 145, 59, 229, 22, 34, 2, 56, 35, 25, 97, 187, 210, 87, 19, …
## $ E_LIMENG   <int> 755, 426, 934, 93, 1076, 36, 664, 175, 129, 447, 461, 509, …
## $ M_LIMENG   <dbl> 209.5, 205.9, 239.3, 137.4, 250.2, 87.9, 215.1, 127.6, 99.3…
## $ E_MUNIT    <int> 67, 886, 211, 134, 1990, 679, 215, 587, 249, 959, 277, 895,…
## $ M_MUNIT    <dbl> 37.1, 308.7, 104.2, 47.4, 303.0, 174.3, 99.9, 147.8, 86.8, …
## $ E_MOBILE   <int> 7770, 4279, 6108, 2625, 7904, 2378, 3140, 2455, 4081, 8615,…
## $ M_MOBILE   <int> 431, 469, 476, 212, 546, 261, 201, 310, 362, 672, 352, 493,…
## $ E_CROWD    <int> 264, 299, 339, 119, 772, 404, 272, 83, 320, 555, 313, 422, …
## $ M_CROWD    <dbl> 77.1, 142.3, 130.7, 57.7, 206.2, 148.9, 86.0, 54.7, 129.1, …
## $ E_NOVEH    <int> 763, 1191, 856, 520, 2599, 989, 1189, 1322, 785, 1279, 1023…
## $ M_NOVEH    <int> 160, 272, 201, 102, 331, 160, 218, 244, 170, 203, 187, 276,…
## $ E_GROUPQ   <int> 654, 546, 543, 322, 3112, 512, 601, 436, 666, 1200, 1263, 4…
## $ M_GROUPQ   <int> 142, 161, 117, 88, 436, 136, 81, 83, 175, 218, 214, 439, 16…
## $ EP_POV     <dbl> -999.0, 15.4, 14.4, 23.5, 18.6, 16.6, 15.1, 16.3, 18.3, 15.…
## $ MP_POV     <dbl> -999.0, 2.1, 1.7, 2.6, 1.2, 2.1, 1.4, 1.6, 1.9, 1.6, 1.8, 1…
## $ EP_UNEMP   <dbl> -999.0, 4.2, 4.1, 6.7, 8.8, 5.0, 5.9, 5.9, 8.7, 5.5, 8.8, 6…
## $ MP_UNEMP   <dbl> -999.0, 1.0, 0.8, 1.7, 1.0, 1.1, 1.0, 1.2, 1.5, 1.0, 1.4, 1…
## $ EP_PCI     <dbl> -999, 29372, 22656, 20430, 24706, 22827, 27577, 24918, 2307…
## $ MP_PCI     <dbl> -999, 2306, 905, 1258, 758, 1482, 976, 1055, 854, 906, 925,…
## $ EP_NOHSDP  <dbl> 13.8, 11.3, 19.8, 15.4, 15.9, 18.6, 13.8, 15.6, 17.1, 18.2,…
## $ MP_NOHSDP  <dbl> 1.6, 1.3, 1.8, 1.9, 1.0, 1.5, 1.2, 1.4, 1.6, 1.4, 1.4, 1.2,…
## $ EP_AGE65   <dbl> 18.0, 14.6, 17.8, 19.0, 16.8, 18.9, 16.3, 19.3, 20.5, 18.0,…
## $ MP_AGE65   <dbl> 0.1, 0.1, 0.2, 0.1, 0.1, 0.3, 0.1, 0.1, 0.3, 0.1, 0.2, 0.1,…
## $ EP_AGE17   <dbl> 23.7, 24.2, 23.4, 22.8, 21.9, 20.7, 23.7, 21.3, 21.9, 22.6,…
## $ MP_AGE17   <dbl> 0.0, 0.1, 0.1, 0.4, 0.1, 0.3, 0.0, 0.1, 0.2, 0.2, 0.1, 0.0,…
## $ EP_DISABL  <dbl> 16.1, 19.3, 14.2, 17.7, 20.8, 16.7, 17.7, 18.7, 20.9, 17.6,…
## $ MP_DISABL  <dbl> 1.3, 1.3, 1.0, 1.9, 0.9, 1.2, 1.0, 1.0, 1.1, 1.2, 1.4, 1.1,…
## $ EP_SNGPNT  <dbl> 10.7, 7.5, 7.0, 10.5, 10.4, 9.7, 9.7, 8.2, 9.0, 6.7, 8.2, 9…
## $ MP_SNGPNT  <dbl> 2.3, 1.5, 1.3, 2.1, 1.0, 1.8, 1.4, 1.3, 1.3, 0.8, 1.3, 1.4,…
## $ EP_MINRTY  <dbl> 87.5, 25.0, 12.9, 48.1, 27.5, 44.2, 29.4, 21.3, 16.6, 7.9, …
## $ MP_MINRTY  <dbl> 0.4, 0.1, 0.4, 0.1, 0.0, 0.0, 0.1, 0.1, 0.1, 0.1, 0.4, 0.3,…
## $ EP_LIMENG  <dbl> 2.1, 0.8, 1.7, 0.5, 1.0, 0.1, 1.4, 0.3, 0.4, 0.6, 1.0, 0.7,…
## $ MP_LIMENG  <dbl> 0.6, 0.4, 0.4, 0.7, 0.2, 0.3, 0.4, 0.2, 0.3, 0.2, 0.5, 0.3,…
## $ EP_MUNIT   <dbl> 0.3, 3.8, 0.9, 1.3, 3.7, 4.0, 0.9, 2.2, 1.3, 2.5, 1.2, 2.6,…
## $ MP_MUNIT   <dbl> 0.2, 1.3, 0.4, 0.5, 0.6, 1.0, 0.4, 0.6, 0.5, 0.5, 0.4, 0.8,…
## $ EP_MOBILE  <dbl> 38.8, 18.4, 25.2, 26.2, 14.7, 14.0, 13.6, 9.3, 21.6, 22.9, …
## $ MP_MOBILE  <dbl> 2.2, 2.0, 2.0, 2.1, 1.0, 1.5, 0.9, 1.2, 1.9, 1.8, 1.5, 1.5,…
## $ EP_CROWD   <dbl> 2.1, 1.4, 1.6, 1.8, 1.7, 3.0, 1.4, 0.4, 2.1, 1.8, 1.7, 1.4,…
## $ MP_CROWD   <dbl> 0.6, 0.7, 0.6, 0.9, 0.5, 1.1, 0.4, 0.3, 0.9, 0.5, 0.6, 0.6,…
## $ EP_NOVEH   <dbl> 6.2, 5.6, 4.2, 7.8, 5.8, 7.3, 6.0, 6.1, 5.2, 4.2, 5.5, 4.7,…
## $ MP_NOVEH   <dbl> 1.3, 1.3, 1.0, 1.5, 0.7, 1.2, 1.1, 1.1, 1.1, 0.7, 1.0, 0.9,…
## $ EP_GROUPQ  <dbl> 1.7, 1.0, 0.9, 1.6, 2.7, 1.5, 1.2, 0.8, 1.8, 1.5, 2.6, 5.4,…
## $ MP_GROUPQ  <dbl> 0.4, 0.3, 0.2, 0.4, 0.4, 0.4, 0.2, 0.2, 0.5, 0.3, 0.4, 0.5,…
## $ EPL_POV    <dbl> -999.0000, 0.5401, 0.4723, 0.8860, 0.7322, 0.6258, 0.5229, …
## $ EPL_UNEMP  <dbl> -999.0000, 0.2745, 0.2611, 0.6968, 0.8850, 0.4166, 0.5669, …
## $ EPL_PCI    <dbl> -999.0000, 0.2860, 0.7561, 0.8879, 0.6076, 0.7465, 0.4022, …
## $ EPL_NOHSDP <dbl> 0.5922, 0.4397, 0.8405, 0.6753, 0.6988, 0.7966, 0.5922, 0.6…
## $ SPL_THEME1 <dbl> -999.0000, 1.5403, 2.3300, 3.1460, 2.9236, 2.5855, 2.0842, …
## $ RPL_THEME1 <dbl> -999.0000, 0.3631, 0.6143, 0.8455, 0.7866, 0.6901, 0.5373, …
## $ EPL_AGE65  <dbl> 0.4893, 0.1850, 0.4715, 0.5928, 0.3664, 0.5817, 0.3200, 0.6…
## $ EPL_AGE17  <dbl> 0.6826, 0.7529, 0.6406, 0.5578, 0.4323, 0.2846, 0.6826, 0.3…
## $ EPL_DISABL <dbl> 0.5610, 0.7905, 0.3763, 0.6845, 0.8564, 0.6074, 0.6845, 0.7…
## $ EPL_SNGPNT <dbl> 0.8383, 0.3792, 0.2961, 0.8185, 0.8077, 0.7316, 0.7316, 0.5…
## $ SPL_THEME2 <dbl> 2.5712, 2.1076, 1.7845, 2.6536, 2.4628, 2.2053, 2.4187, 2.2…
## $ RPL_THEME2 <dbl> 0.8758, 0.5810, 0.3187, 0.9077, 0.8303, 0.6609, 0.8067, 0.6…
## $ EPL_MINRTY <dbl> 0.9917, 0.6336, 0.4206, 0.8711, 0.6616, 0.8437, 0.6883, 0.5…
## $ EPL_LIMENG <dbl> 0.7740, 0.5113, 0.7170, 0.3582, 0.5791, 0.0707, 0.6718, 0.2…
## $ SPL_THEME3 <dbl> 1.7657, 1.1449, 1.1376, 1.2293, 1.2407, 0.9144, 1.3601, 0.8…
## $ RPL_THEME3 <dbl> 0.9268, 0.5947, 0.5915, 0.6447, 0.6507, 0.4597, 0.7256, 0.3…
## $ EPL_MUNIT  <dbl> 0.0551, 0.6017, 0.1512, 0.2416, 0.5931, 0.6240, 0.1512, 0.4…
## $ EPL_MOBILE <dbl> 0.9869, 0.7408, 0.8816, 0.8940, 0.6396, 0.6192, 0.6065, 0.4…
## $ EPL_CROWD  <dbl> 0.5498, 0.2964, 0.3703, 0.4457, 0.4056, 0.7622, 0.2964, 0.0…
## $ EPL_NOVEH  <dbl> 0.5788, 0.4846, 0.2420, 0.7685, 0.5142, 0.7224, 0.5441, 0.5…
## $ EPL_GROUPQ <dbl> 0.4126, 0.1525, 0.1165, 0.3792, 0.6250, 0.3457, 0.2184, 0.0…
## $ SPL_THEME4 <dbl> 2.5832, 2.2760, 1.7616, 2.7290, 2.7775, 3.0735, 1.8166, 1.5…
## $ RPL_THEME4 <dbl> 0.5409, 0.3741, 0.1741, 0.6259, 0.6492, 0.7943, 0.1923, 0.1…
## $ SPL_THEMES <dbl> -999.0000, 7.0688, 7.0137, 9.7579, 9.4046, 8.7787, 7.6796, …
## $ RPL_THEMES <dbl> -999.0000, 0.4354, 0.4242, 0.8653, 0.8252, 0.7382, 0.5408, …
## $ F_POV      <int> -999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_UNEMP    <int> -999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_PCI      <int> -999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_NOHSDP   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_THEME1   <int> -999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_AGE65    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_AGE17    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_DISABL   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_SNGPNT   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_THEME2   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_MINRTY   <int> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_LIMENG   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_THEME3   <int> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_MUNIT    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_MOBILE   <int> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_CROWD    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_NOVEH    <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_GROUPQ   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_THEME4   <int> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ F_TOTAL    <int> -999, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ E_UNINSUR  <int> 4160, 3875, 6303, 2005, 10686, 3613, 5004, 4353, 4116, 9631…
## $ M_UNINSUR  <int> 588, 508, 732, 340, 796, 476, 540, 552, 506, 933, 600, 816,…
## $ EP_UNINSUR <dbl> 10.6, 7.1, 11.0, 10.2, 9.4, 10.8, 10.1, 8.1, 11.2, 11.8, 11…
## $ MP_UNINSUR <dbl> 1.5, 0.9, 1.3, 1.7, 0.7, 1.4, 1.1, 1.0, 1.4, 1.1, 1.3, 1.1,…
## $ E_DAYPOP   <int> 32290, 37301, 40036, 17280, 117894, 27176, 44908, 56580, 35…
dim(svi_2018)
## [1] 3142  123

The USGS glyphosate data has 15314 rows and 6 columns. The US County shapefile data has 3221 rows and 10 columns. The PRISM county-level data for 2013 has 1134420 rows and 6 columns. These data are connected by FIPS code, which is comprised of a two-digit state code and then a three digit county code. The USGS pesticide estimate data contains a low and high estimate (EPEST_LOW_KG and EPEST_HIGH_KG). The CDC SVI data has 3142 rows and 133 columns.

For the SVI, the four summary theme ranking variables, detailed in the Data Dictionary below, are:

Socioeconomic theme – RPL_THEME1 Household Composition and Disability – RPL_THEME2 Minority Status & Language – RPL_THEME3 Housing & Transportation – RPL_THEME4 The overall tract summary ranking variable is RPL_THEMES.

Variables that start with E are the estimates and M are the margins of error. Full documentation for SVI data can be found here: https://www.atsdr.cdc.gov/placeandhealth/svi/documentation/SVI_documentation_2014.html.

Data Cleaning/Wrangling

Create combined FIPS code (two-digit state plus three-digit county in format 01001).

usgs_gly_2013 <- usgs_gly_2013 %>% mutate(STATE_FIPS_CODE =
                              str_pad(as.character(usgs_gly_2013$STATE_FIPS_CODE),
                              2, pad = "0"),
                COUNTY_FIPS_CODE =                  
                  str_pad(as.character(usgs_gly_2013$COUNTY_FIPS_CODE),
                              3, pad = "0"),
                fips = as.factor(paste0(STATE_FIPS_CODE, COUNTY_FIPS_CODE)))

usgs_gly_2014 <- usgs_gly_2014 %>% mutate(STATE_FIPS_CODE =
                              str_pad(as.character(usgs_gly_2014$STATE_FIPS_CODE),
                              2, pad = "0"),
                COUNTY_FIPS_CODE =                  
                  str_pad(as.character(usgs_gly_2014$COUNTY_FIPS_CODE),
                              3, pad = "0"),
                fips = as.factor(paste0(STATE_FIPS_CODE, COUNTY_FIPS_CODE)))

usgs_gly_2015 <- usgs_gly_2015 %>% mutate(STATE_FIPS_CODE =
                              str_pad(as.character(usgs_gly_2015$STATE_FIPS_CODE),
                              2, pad = "0"),
                COUNTY_FIPS_CODE =                  
                  str_pad(as.character(usgs_gly_2015$COUNTY_FIPS_CODE),
                              3, pad = "0"),
                fips = as.factor(paste0(STATE_FIPS_CODE, COUNTY_FIPS_CODE)))

usgs_gly_2016 <- usgs_gly_2016 %>% mutate(STATE_FIPS_CODE =
                              str_pad(as.character(usgs_gly_2016$STATE_FIPS_CODE),
                              2, pad = "0"),
                COUNTY_FIPS_CODE =                  
                  str_pad(as.character(usgs_gly_2016$COUNTY_FIPS_CODE),
                              3, pad = "0"),
                fips = as.factor(paste0(STATE_FIPS_CODE, COUNTY_FIPS_CODE)))

usgs_gly_2017 <- usgs_gly_2017 %>% mutate(STATE_FIPS_CODE =
                              str_pad(as.character(usgs_gly_2017$STATE_FIPS_CODE),
                              2, pad = "0"),
                COUNTY_FIPS_CODE =                  
                  str_pad(as.character(usgs_gly_2017$COUNTY_FIPS_CODE),
                              3, pad = "0"),
                fips = as.factor(paste0(STATE_FIPS_CODE, COUNTY_FIPS_CODE)))

us_county_2013 <- us_county_2013 %>% mutate(fips = 
                                  as.factor(paste0(STATEFP, COUNTYFP)))

svi_2014 <- svi_2014 %>% mutate(fips =
                              str_pad(as.character(svi_2014$FIPS),
                              5, pad = "0"))

svi_2016 <- svi_2016 %>% mutate(fips =
                              str_pad(as.character(svi_2016$FIPS),
                              5, pad = "0"))

svi_2018 <- svi_2018 %>% mutate(fips =
                              str_pad(as.character(svi_2018$FIPS),
                              5, pad = "0"))
#2013
prism_county_2013 <- prism_county_2013 %>% mutate(date=as.Date(date,
                                              format="%d/%m/%Y")) %>%
  select(-c(day,month, year))

group_by(prism_county_2013, fips) %>% tally()
## # A tibble: 3,108 × 2
##    fips      n
##    <fct> <int>
##  1 01001   365
##  2 01003   365
##  3 01005   365
##  4 01007   365
##  5 01009   365
##  6 01011   365
##  7 01013   365
##  8 01015   365
##  9 01017   365
## 10 01019   365
## # ℹ 3,098 more rows
#got 3108 fips codes

prism_01001 <- prism_county_2013 %>% filter(fips=="01001")
quantile(prism_01001$tmean, probs = 0.85)
##      85% 
## 26.03194
prism_12033 <- prism_county_2013 %>% filter(fips=="12033")
quantile(prism_12033$tmean, probs = 0.85)
##      85% 
## 26.83704
prism_county_2013_85s <- prism_county_2013 %>% group_by(fips) %>% mutate(tmean85_fullyear = quantile(tmean, probs=0.85, na.rm=TRUE)) %>% 
  filter(between(date, as.Date('2013-04-30'), as.Date('2013-09-01'))) %>%
  mutate(tmean85_spray = quantile(tmean, probs=0.85, na.rm=TRUE)) %>% 
  summarise(tmean85_fullyear=mean(tmean85_fullyear),
            tmean85_spray=mean(tmean85_spray))

#2014
prism_county_2014 <- prism_county_2014 %>% mutate(date=as.Date(date,
                                              format="%d/%m/%Y")) %>%
  select(-c(day,month, year))

group_by(prism_county_2014, fips) %>% tally()
## # A tibble: 3,108 × 2
##    fips      n
##    <fct> <int>
##  1 01001   365
##  2 01003   365
##  3 01005   365
##  4 01007   365
##  5 01009   365
##  6 01011   365
##  7 01013   365
##  8 01015   365
##  9 01017   365
## 10 01019   365
## # ℹ 3,098 more rows
#got 3108 fips codes

prism_01001 <- prism_county_2014 %>% filter(fips=="01001")
quantile(prism_01001$tmean, probs = 0.85)
##      85% 
## 26.47876
prism_12033 <- prism_county_2014 %>% filter(fips=="12033")
quantile(prism_12033$tmean, probs = 0.85)
##      85% 
## 27.25865
prism_county_2014_85s <- prism_county_2014 %>% group_by(fips) %>% mutate(tmean85_fullyear = quantile(tmean, probs=0.85, na.rm=TRUE)) %>% 
  filter(between(date, as.Date('2014-04-30'), as.Date('2014-09-01'))) %>%
  mutate(tmean85_spray = quantile(tmean, probs=0.85, na.rm=TRUE)) %>% 
  summarise(tmean85_fullyear=mean(tmean85_fullyear),
            tmean85_spray=mean(tmean85_spray))

#2015
prism_county_2015 <- prism_county_2015 %>% mutate(date=as.Date(date,
                                              format="%d/%m/%Y")) %>%
  select(-c(day,month, year))

group_by(prism_county_2015, fips) %>% tally()
## # A tibble: 3,108 × 2
##    fips      n
##    <fct> <int>
##  1 01001   365
##  2 01003   365
##  3 01005   365
##  4 01007   365
##  5 01009   365
##  6 01011   365
##  7 01013   365
##  8 01015   365
##  9 01017   365
## 10 01019   365
## # ℹ 3,098 more rows
#got 3108 fips codes

prism_county_2015_85s <- prism_county_2015 %>% group_by(fips) %>% mutate(tmean85_fullyear = quantile(tmean, probs=0.85, na.rm=TRUE)) %>% 
  filter(between(date, as.Date('2015-04-30'), as.Date('2015-09-01'))) %>%
  mutate(tmean85_spray = quantile(tmean, probs=0.85, na.rm=TRUE)) %>% 
  summarise(tmean85_fullyear=mean(tmean85_fullyear),
            tmean85_spray=mean(tmean85_spray))

#2016
prism_county_2016 <- prism_county_2016 %>% mutate(date=as.Date(date,
                                              format="%d/%m/%Y")) %>%
  select(-c(day,month, year))

group_by(prism_county_2016, fips) %>% tally()
## # A tibble: 3,108 × 2
##    fips      n
##    <fct> <int>
##  1 01001   366
##  2 01003   366
##  3 01005   366
##  4 01007   366
##  5 01009   366
##  6 01011   366
##  7 01013   366
##  8 01015   366
##  9 01017   366
## 10 01019   366
## # ℹ 3,098 more rows
#got 3108 fips codes

prism_county_2016_85s <- prism_county_2016 %>% group_by(fips) %>% mutate(tmean85_fullyear = quantile(tmean, probs=0.85, na.rm=TRUE)) %>% 
  filter(between(date, as.Date('2016-04-30'), as.Date('2016-09-01'))) %>%
  mutate(tmean85_spray = quantile(tmean, probs=0.85, na.rm=TRUE)) %>% 
  summarise(tmean85_fullyear=mean(tmean85_fullyear),
            tmean85_spray=mean(tmean85_spray))

#2017
prism_county_2017 <- prism_county_2017 %>% mutate(date=as.Date(date,
                                              format="%d/%m/%Y")) %>%
  select(-c(day,month, year))

group_by(prism_county_2017, fips) %>% tally()
## # A tibble: 3,108 × 2
##    fips      n
##    <fct> <int>
##  1 01001   365
##  2 01003   365
##  3 01005   365
##  4 01007   365
##  5 01009   365
##  6 01011   365
##  7 01013   365
##  8 01015   365
##  9 01017   365
## 10 01019   365
## # ℹ 3,098 more rows
#got 3108 fips codes

prism_county_2017_85s <- prism_county_2017 %>% group_by(fips) %>% mutate(tmean85_fullyear = quantile(tmean, probs=0.85, na.rm=TRUE)) %>% 
  filter(between(date, as.Date('2017-04-30'), as.Date('2017-09-01'))) %>%
  mutate(tmean85_spray = quantile(tmean, probs=0.85, na.rm=TRUE)) %>% 
  summarise(tmean85_fullyear=mean(tmean85_fullyear),
            tmean85_spray=mean(tmean85_spray))
svi_2014 <- svi_2014 %>% select(c(AFFGEOID, ST, STATE, ST_ABBR, COUNTY,fips,
                                  RPL_THEME1, RPL_THEME2, RPL_THEME3,
                                  RPL_THEME4, RPL_THEMES, Shape,
                                  "Shape.STArea()","Shape.STLength()")) %>% 
  rename(SVI_socioecon=RPL_THEME1, SVI_housecompdis=RPL_THEME2, 
         SVI_minstatlang=RPL_THEME3, SVI_housing=RPL_THEME4,
         SVI_overall=RPL_THEMES)

svi_2016 <- svi_2016 %>% select(c(ST, STATE, ST_ABBR, COUNTY,fips,
                                  RPL_THEME1, RPL_THEME2, RPL_THEME3,
                                  RPL_THEME4, RPL_THEMES)) %>% 
  rename(SVI_socioecon=RPL_THEME1, SVI_housecompdis=RPL_THEME2, 
         SVI_minstatlang=RPL_THEME3, SVI_housing=RPL_THEME4,
         SVI_overall=RPL_THEMES)

#found a -999 value in svi_2018, filtering out
svi_2018 <- svi_2018 %>% select(c(ST, STATE, ST_ABBR, COUNTY,fips,
                                  RPL_THEME1, RPL_THEME2, RPL_THEME3,
                                  RPL_THEME4, RPL_THEMES)) %>% 
  rename(SVI_socioecon=RPL_THEME1, SVI_housecompdis=RPL_THEME2, 
         SVI_minstatlang=RPL_THEME3, SVI_housing=RPL_THEME4,
         SVI_overall=RPL_THEMES) %>% filter(SVI_overall!=-999)

svi_2014 %>% summarize(mean_svi=mean(SVI_overall), sd_svi=sd(SVI_overall),
                       max_svi=max(SVI_overall), min_svi=min(SVI_overall),
                       n=length(unique(fips)))
##    mean_svi    sd_svi max_svi min_svi    n
## 1 0.4999994 0.2888131       1       0 3142
svi_2016 %>% summarize(mean_svi=mean(SVI_overall), sd_svi=sd(SVI_overall),
                       max_svi=max(SVI_overall), min_svi=min(SVI_overall),
                       n=length(unique(fips)))
##    mean_svi   sd_svi max_svi min_svi    n
## 1 0.4999993 0.288813       1       0 3142
svi_2018 %>% summarize(mean_svi=mean(SVI_overall), sd_svi=sd(SVI_overall),
                       max_svi=max(SVI_overall), min_svi=min(SVI_overall),
                       n=length(unique(fips)))
##   mean_svi    sd_svi max_svi min_svi    n
## 1 0.499993 0.2888138       1       0 3141
merged_2013 <- left_join(us_county_2013, prism_county_2013_85s) %>% 
  left_join(usgs_gly_2013) %>% filter(STATEFP != "02" & STATEFP != "15" &
                                   STATEFP != "72") %>%
  mutate(gly_per_land_low = EPEST_LOW_KG/(ALAND/1e+6),
         gly_per_land_high = EPEST_HIGH_KG/(ALAND/1e+6),
         gly_per_totalarea_low = EPEST_LOW_KG/((ALAND+AWATER)/1e+6),
         gly_per_totalarea_high = EPEST_HIGH_KG/((ALAND+AWATER)/1e+6)) %>%
  left_join(svi_2014)
## Joining with `by = join_by(fips)`
## Joining with `by = join_by(fips)`
## Joining with `by = join_by(AFFGEOID, fips)`
merged_2014 <- left_join(us_county_2013, prism_county_2014_85s) %>% 
  left_join(usgs_gly_2014) %>% filter(STATEFP != "02" & STATEFP != "15" &
                                   STATEFP != "72") %>%
  mutate(gly_per_land_low = EPEST_LOW_KG/(ALAND/1e+6),
         gly_per_land_high = EPEST_HIGH_KG/(ALAND/1e+6),
         gly_per_totalarea_low = EPEST_LOW_KG/((ALAND+AWATER)/1e+6),
         gly_per_totalarea_high = EPEST_HIGH_KG/((ALAND+AWATER)/1e+6)) %>%
  left_join(svi_2014)
## Joining with `by = join_by(fips)`
## Joining with `by = join_by(fips)`
## Joining with `by = join_by(AFFGEOID, fips)`
merged_2015 <- left_join(us_county_2013, prism_county_2015_85s) %>% 
  left_join(usgs_gly_2015) %>% filter(STATEFP != "02" & STATEFP != "15" &
                                   STATEFP != "72") %>%
  mutate(gly_per_land_low = EPEST_LOW_KG/(ALAND/1e+6),
         gly_per_land_high = EPEST_HIGH_KG/(ALAND/1e+6),
         gly_per_totalarea_low = EPEST_LOW_KG/((ALAND+AWATER)/1e+6),
         gly_per_totalarea_high = EPEST_HIGH_KG/((ALAND+AWATER)/1e+6)) %>%
  left_join(svi_2016)
## Joining with `by = join_by(fips)`
## Joining with `by = join_by(fips)`
## Joining with `by = join_by(fips)`
merged_2016 <- left_join(us_county_2013, prism_county_2016_85s) %>% 
  left_join(usgs_gly_2016) %>% filter(STATEFP != "02" & STATEFP != "15" &
                                   STATEFP != "72") %>%
  mutate(gly_per_land_low = EPEST_LOW_KG/(ALAND/1e+6),
         gly_per_land_high = EPEST_HIGH_KG/(ALAND/1e+6),
         gly_per_totalarea_low = EPEST_LOW_KG/((ALAND+AWATER)/1e+6),
         gly_per_totalarea_high = EPEST_HIGH_KG/((ALAND+AWATER)/1e+6)) %>%
  left_join(svi_2016)
## Joining with `by = join_by(fips)`
## Joining with `by = join_by(fips)`
## Joining with `by = join_by(fips)`
merged_2017 <- left_join(us_county_2013, prism_county_2017_85s) %>% 
  left_join(usgs_gly_2017) %>% filter(STATEFP != "02" & STATEFP != "15" &
                                   STATEFP != "72") %>%
  mutate(gly_per_land_low = EPEST_LOW_KG/(ALAND/1e+6),
         gly_per_land_high = EPEST_HIGH_KG/(ALAND/1e+6),
         gly_per_totalarea_low = EPEST_LOW_KG/((ALAND+AWATER)/1e+6),
         gly_per_totalarea_high = EPEST_HIGH_KG/((ALAND+AWATER)/1e+6)) %>%
  left_join(svi_2018)
## Joining with `by = join_by(fips)`
## Joining with `by = join_by(fips)`
## Joining with `by = join_by(fips)`
#cleaned merged data so that it's just for conterminous U.S. (not AK which is 02, HI which is 15, or PR which is )
compounddata_formap <- function(compound, year) {
  usgs_data <- usgs_ests %>% filter(COMPOUND==compound) %>% 
    filter(YEAR==year) %>% mutate(STATE_FIPS_CODE = 
                        str_pad(as.character(STATE_FIPS_CODE),
                              2, pad = "0"),
                COUNTY_FIPS_CODE =                  
                  str_pad(as.character(COUNTY_FIPS_CODE),
                              3, pad = "0"),
                fips = as.factor(paste0(STATE_FIPS_CODE, 
                                        COUNTY_FIPS_CODE)))
  
  if (year==2013) {
    us_county_data <- us_county_2013
    prism_county_data <- prism_county_2013_85s
    svi_data <- svi_2014
  }
  
  if (year==2014) {
    prism_county_data <- prism_county_2014_85s
    svi_data <- svi_2014
  }
  
  if (year==2015) {
    prism_county_data <- prism_county_2015_85s
    svi_data <- svi_2016
  }
  
  if (year==2016) {
    prism_county_data <- prism_county_2016_85s
    svi_data <- svi_2016
  }
  
  if (year==2017) {
    prism_county_data <- prism_county_2017_85s
    svi_data <- svi_2018
  }
  
  merged_data <- left_join(us_county_2013, prism_county_data) %>% 
  left_join(usgs_data) %>% filter(STATEFP != "02" & STATEFP != "15" &
                                   STATEFP != "72") %>%
  mutate(epest_per_land_low = EPEST_LOW_KG/(ALAND/1e+6),
         epest_per_land_high = EPEST_HIGH_KG/(ALAND/1e+6),
         epest_per_totalarea_low = EPEST_LOW_KG/((ALAND+AWATER)/1e+6),
         epest_per_totalarea_high = EPEST_HIGH_KG/((ALAND+AWATER)/1e+6)) %>%
  left_join(svi_data)

  merged_data
}

hot_compounddata_formap <- function(compound, year, temp_c) {
  compounddata_formap(compound, year) %>% filter(tmean85_spray>=temp_c)
}

Data Visualization

Glyphosate Usage 2013-2017

plot_gly <- function(year, est_type) {
  if (year==2013) {
    data = merged_2013
  }
  if (year==2014) {
    data = merged_2014
   }
  if (year==2015) {
    data = merged_2015
  }
  if (year==2016) {
    data = merged_2016
  }
  if (year==2017) {
    data = merged_2017
  }
  
  if (est_type=="low") {
    map <- ggplot(data) +
    geom_sf(aes(fill = gly_per_land_low), linetype=0) +
    ggtitle("Glyphosate per Land Area (Low Est.)", subtitle=year) +
    theme_void(base_size = 14) +
    theme(
      plot.title = element_text(hjust = 0.5),
      axis.text.x = element_blank(),
      axis.text.y = element_blank(),
      axis.ticks = element_blank(),
      rect = element_blank()
    ) +
    scale_fill_gradientn("Glyphosate usage \n(kg/km^2)",
                         colors = met.brewer("Tam"),
                         na.value = "grey50")
  }
  
  if (est_type=="high") {
    map <- ggplot(data) +
    geom_sf(aes(fill = gly_per_land_high), linetype=0) +
    ggtitle("Glyphosate per Land Area (High Est.)", subtitle=year) +
    theme_void(base_size = 14) +
    theme(
      plot.title = element_text(hjust = 0.5),
      axis.text.x = element_blank(),
      axis.text.y = element_blank(),
      axis.ticks = element_blank(),
      rect = element_blank()
    ) +
    scale_fill_gradientn("Glyphosate usage \n(kg/km^2)",
                         colors = met.brewer("Tam"),
                         na.value = "grey50")
  }
  map
}

plot_gly(2013, "low")
plot_gly(2013, "high")
plot_gly(2014, "low")
plot_gly(2014, "high")
plot_gly(2015, "low")
plot_gly(2015, "high")
plot_gly(2016, "low")
plot_gly(2016, "high")
plot_gly(2017, "low")
plot_gly(2017, "high")

85th Percentile Mean Daily Temperatures 2013-2017

I created plots for temperature during spray season and during the full year.

2013

temp85_plot <- ggplot() +
  geom_sf(data = merged_2013, aes(fill = tmean85_spray), linetype=0) +
  ggtitle("85th Percentile Mean Temp. by County (2013)", 
          subtitle="Conterminous U.S.") +  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Temperature (C)",
                       colors = rev(met.brewer("Hiroshige")),
                       na.value = "grey50")

temp85_plot_f_2013 <- ggplot() +
  geom_sf(data = merged_2013, aes(fill = (tmean85_spray*9/5+32)), linetype=0) +
  ggtitle("85th Percentile Mean Temp. by County (2013)", 
          subtitle="Conterminous U.S.") +
  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Temperature (F)",
                       colors = rev(met.brewer("Hiroshige")),
                       na.value = "grey50")
plot(temp85_plot)
plot(temp85_plot_f_2013)

2014

temp85_plot <- ggplot() +
  geom_sf(data = merged_2014, aes(fill = tmean85_spray), linetype=0) +
  ggtitle("85th Percentile Mean Temp. by County (2014)", 
          subtitle="Conterminous U.S.") +  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Temperature (C)",
                       colors = rev(met.brewer("Hiroshige")),
                       na.value = "grey50")

temp85_plot_f <- ggplot() +
  geom_sf(data = merged_2014, aes(fill = (tmean85_spray*9/5+32)), linetype=0) +
  ggtitle("85th Percentile Mean Temp. by County (2014)", 
          subtitle="Conterminous U.S.") +
  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Temperature (F)",
                       colors = rev(met.brewer("Hiroshige")),
                       na.value = "grey50")
plot(temp85_plot)
plot(temp85_plot_f)

plot(temp85_plot)
plot_gly(2014, "high")

2015

temp85_plot <- ggplot() +
  geom_sf(data = merged_2015, aes(fill = tmean85_spray), linetype=0) +
  ggtitle("85th Percentile Mean Temp. by County (2015)", 
          subtitle="Conterminous U.S.") +  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Temperature (C)",
                       colors = rev(met.brewer("Hiroshige")),
                       na.value = "grey50")

temp85_plot_f <- ggplot() +
  geom_sf(data = merged_2015, aes(fill = (tmean85_spray*9/5+32)), linetype=0) +
  ggtitle("85th Percentile Mean Temp. by County (2015)", 
          subtitle="Conterminous U.S.") +
  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Temperature (F)",
                       colors = rev(met.brewer("Hiroshige")),
                       na.value = "grey50")
plot(temp85_plot)
plot(temp85_plot_f)

2016

temp85_plot <- ggplot() +
  geom_sf(data = merged_2016, aes(fill = tmean85_spray), linetype=0) +
  ggtitle("85th Percentile Mean Temp. by County (2016)", 
          subtitle="Conterminous U.S.") +  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Temperature (C)",
                       colors = rev(met.brewer("Hiroshige")),
                       na.value = "grey50")

temp85_plot_f <- ggplot() +
  geom_sf(data = merged_2016, aes(fill = (tmean85_spray*9/5+32)), linetype=0) +
  ggtitle("85th Percentile Mean Temp. by County (2016)", 
          subtitle="Conterminous U.S.") +
  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Temperature (F)",
                       colors = rev(met.brewer("Hiroshige")),
                       na.value = "grey50")
plot(temp85_plot)
plot(temp85_plot_f)

2017

temp85_plot <- ggplot() +
  geom_sf(data = merged_2017, aes(fill = tmean85_spray), linetype=0) +
  ggtitle("85th Percentile Mean Temp. by County (2017)", 
          subtitle="Conterminous U.S.") +  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Temperature (C)",
                       colors = rev(met.brewer("Hiroshige")),
                       na.value = "grey50")

temp85_plot_f_2017 <- ggplot() +
  geom_sf(data = merged_2017, aes(fill = (tmean85_spray*9/5+32)), linetype=0) +
  ggtitle("85th Percentile Mean Temp. by County (2017)", 
          subtitle="Conterminous U.S.") +
  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Temperature (F)",
                       colors = rev(met.brewer("Hiroshige")),
                       na.value = "grey50")
plot(temp85_plot)
plot(temp85_plot_f)

plot(temp85_plot)
plot_gly(2017, "high")

Higher SVI Counties

# Function to create SVI map
map_svi <- function(data, svi_var, year, percentile) {
  # Filter data for the given year
  filtered_data <- data %>%
    filter(YEAR == year) %>%
    filter(!!sym(svi_var) >= percentile)  # Use rlang's !!sym to evaluate svi_var correctly
  
  # Create the map
  ggplot() +
    # Layer for counties with high percentile (bold black outline)
    geom_sf(
      data = filtered_data,
      aes(fill = !!sym(svi_var)),  # Ensure the svi_var is evaluated correctly in aes
      size = 2.0,  # Thicker border
      color = "black"
    ) +
    # Layer for counties below the percentile (no fill or outline)
    geom_sf(
      data = data %>%
        filter(YEAR == year) %>%
        filter(!!sym(svi_var) < percentile),
      aes(),
      fill = NA,
      color = NA
    ) +
    scale_fill_viridis_c() +  # Using viridis color scale for fill
    labs(
      title = paste("Map of", svi_var, "in", year),
      fill = svi_var
    ) +
    theme_minimal()
}
# Example of usage
# create_svi_map(data = my_data, svi_var = "SVI", year = 2023, percentile = 90, region_shapefile = "path/to/shapefile.shp")

#housing map for 2017
map_svi(merged_2017, "SVI_housing", 2017, 0.749)

#socioeconomic SVI 2017
map_svi(merged_2017, "SVI_socioecon", 2017, 0.749)

#household composition and disability 2017
map_svi(merged_2017,"SVI_housecompdis", 2017, 0.749)

#minority status and language 2017
map_svi(merged_2017, "SVI_minstatlang", 2017, 0.749)

#overall
map_svi(merged_2017, "SVI_overall", 2017, 0.749)

Data Analysis

I found which counties are hottest during the spray season (85th percentile of mean daily temperature 80 degrees F or hotter), and mapped the glyphosate usage for these counties.

Glyphosate Usage in Hot Counties 2013

all_leaflet_low <- function(compound, year) {

  #get data
  data <- compounddata_formap(compound, year)
  
  map <- data %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Voyager") %>%
    addPolygons(
             stroke = FALSE,
             smoothFactor = 0,
             fillOpacity = 0.5,
             fillColor= ~
               colorNumeric("YlOrRd", epest_per_land_low)(epest_per_land_low)) %>%
  addLegend(pal = colorNumeric("YlOrRd", data$epest_per_land_low),  # Define color palette
        values = data$epest_per_land_low,  # Values to map
        title = "Pesticide Usage (kg/km^2)",  # Legend title
        position = "bottomright",  # Position of the legend
        na.label = "No data"
    )
  
  map
}

all_leaflet_high <- function(compound, year) {

  #get data
  data <- compounddata_formap(compound, year)
  
  map <- data %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Voyager") %>%
    addPolygons(
             stroke = FALSE,
             smoothFactor = 0,
             fillOpacity = 0.5,
             fillColor= ~
               colorNumeric("YlOrRd", epest_per_land_high)(epest_per_land_high)) %>%
  addLegend(pal = colorNumeric("YlOrRd", data$epest_per_land_high),  # Define color palette
        values = data$epest_per_land_high,  # Values to map
        title = "Pesticide Usage (kg/km^2)",  # Legend title
        position = "bottomright",  # Position of the legend
        na.label = "No data"
    )
  
  map
}

hot_leaflet_low <- function(compound, year, temp_c) {

  #get data
  data <- hot_compounddata_formap(compound, year, temp_c)
  
  map <- data %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Voyager") %>%
    addPolygons(
             stroke = FALSE,
             smoothFactor = 0,
             fillOpacity = 0.5,
             fillColor= ~
               colorNumeric("YlOrRd", epest_per_land_low)(epest_per_land_low)) %>%
  addLegend(pal = colorNumeric("YlOrRd", data$epest_per_land_low),  # Define color palette
        values = data$epest_per_land_low,  # Values to map
        title = "Pesticide Usage (kg/km^2)",  # Legend title
        position = "bottomright",  # Position of the legend
        na.label = "No data"
    )
  
  map
}

hot_leaflet_high <- function(compound, year, temp_c) {

  #get data
  data <- hot_compounddata_formap(compound, year, temp_c)
  
  map <- data %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Voyager") %>%
    addPolygons(
             stroke = FALSE,
             smoothFactor = 0,
             fillOpacity = 0.5,
             fillColor= ~
               colorNumeric("YlOrRd", epest_per_land_high)(epest_per_land_high)) %>%
  addLegend(pal = colorNumeric("YlOrRd", data$epest_per_land_high),  # Define color palette
        values = data$epest_per_land_high,  # Values to map
        title = "Pesticide Usage (kg/km^2)",  # Legend title
        position = "bottomright",  # Position of the legend
        na.label = "No data"
    )
  
  map
}
#find areas that have 85th percentile mean daily temp of 80 degrees F or higher
hot_merged_2013 <- merged_2013 %>% filter(tmean85_spray>=26.66)
hot_plot <- ggplot() +
  geom_sf(data = hot_merged_2013, aes(fill = gly_per_land_low), linetype=0) +
  ggtitle("Glyphosate per Land Area", 
          subtitle="Hot Counties") +  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Glyphosate usage",
                       colors = met.brewer("Tam"),
                       na.value = "grey50")
plot(hot_plot)

hot_leaflet <- hot_merged_2013 %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Voyager") %>%
    addPolygons(
             stroke = FALSE,
             smoothFactor = 0,
             fillOpacity = 0.5,
             fillColor= ~
               colorNumeric("YlOrRd",gly_per_land_low)(gly_per_land_low)) %>%
  addLegend(pal = colorNumeric("YlOrRd", hot_merged_2013$gly_per_land_low),  # Define color palette
        values = hot_merged_2013$gly_per_land_low,  # Values to map
        title = "Glyphosate Usage (kg/km^2)",  # Legend title
        position = "bottomright",  # Position of the legend
        na.label = "No data"
    )
## Warning in CPL_crs_from_input(x): GDAL Message 1: +init=epsg:XXXX syntax is
## deprecated. It might return a CRS with a non-EPSG compliant axis order.
hot_leaflet

Glyphosate Usage in Hot Counties 2014

#find areas that have 85th percentile mean daily temp of 80 degrees F or higher
hot_merged_2014 <- merged_2014 %>% filter(tmean85_spray>=26.66)
hot_plot <- ggplot() +
  geom_sf(data = hot_merged_2014, aes(fill = gly_per_land_low), linetype=0) +
  ggtitle("Glyphosate per Land Area", 
          subtitle="Hot Counties") +  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Glyphosate usage",
                       colors = met.brewer("Tam"),
                       na.value = "grey50")
plot(hot_plot)

hot_leaflet <- hot_merged_2014 %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Voyager") %>%
    addPolygons(
             stroke = FALSE,
             smoothFactor = 0,
             fillOpacity = 0.5,
             fillColor= ~
               colorNumeric("YlOrRd",gly_per_land_low)(gly_per_land_low)) %>%
  addLegend(pal = colorNumeric("YlOrRd", hot_merged_2014$gly_per_land_low),  # Define color palette
        values = hot_merged_2014$gly_per_land_low,  # Values to map
        title = "Glyphosate Usage (kg/km^2)",  # Legend title
        position = "bottomright"  # Position of the legend
    )
hot_leaflet

Glyphosate Usage in Hot Counties 2015

#find areas that have 85th percentile mean daily temp of 80 degrees F or higher
hot_merged_2015 <- merged_2015 %>% filter(tmean85_spray>=26.66)
hot_plot <- ggplot() +
  geom_sf(data = hot_merged_2015, aes(fill = gly_per_land_low), linetype=0) +
  ggtitle("Glyphosate per Land Area", 
          subtitle="Hot Counties") +  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Glyphosate usage",
                       colors = met.brewer("Tam"),
                       na.value = "grey50")
plot(hot_plot)

hot_leaflet <- hot_merged_2015 %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Voyager") %>%
    addPolygons(
             stroke = FALSE,
             smoothFactor = 0,
             fillOpacity = 0.5,
             fillColor= ~
               colorNumeric("YlOrRd",gly_per_land_low)(gly_per_land_low)) %>%
  addLegend(pal = colorNumeric("YlOrRd", hot_merged_2015$gly_per_land_low),  # Define color palette
        values = hot_merged_2015$gly_per_land_low,  # Values to map
        title = "Glyphosate Usage (kg/km^2)",  # Legend title
        position = "bottomright"  # Position of the legend
    )
hot_leaflet

Glyphosate Usage in Hot Counties 2016

#find areas that have 85th percentile mean daily temp of 80 degrees F or higher
hot_merged_2016 <- merged_2016 %>% filter(tmean85_spray>=26.66)
hot_plot <- ggplot() +
  geom_sf(data = hot_merged_2016, aes(fill = gly_per_land_low), linetype=0) +
  ggtitle("Glyphosate per Land Area", 
          subtitle="Hot Counties") +  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Glyphosate usage",
                       colors = met.brewer("Tam"),
                       na.value = "grey50")
plot(hot_plot)

hot_leaflet <- hot_merged_2016 %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Voyager") %>%
    addPolygons(
             stroke = FALSE,
             smoothFactor = 0,
             fillOpacity = 0.5,
             fillColor= ~
               colorNumeric("YlOrRd",gly_per_land_low)(gly_per_land_low)) %>%
  addLegend(pal = colorNumeric("YlOrRd", hot_merged_2016$gly_per_land_low),  # Define color palette
        values = hot_merged_2016$gly_per_land_low,  # Values to map
        title = "Glyphosate Usage (kg/km^2)",  # Legend title
        position = "bottomright"  # Position of the legend
    )
hot_leaflet

Glyphosate Usage in Hot Counties 2017

#find areas that have 85th percentile mean daily temp of 80 degrees F or higher
hot_merged_2017 <- merged_2017 %>% filter(tmean85_spray>=26.66)
hot_plot <- ggplot() +
  geom_sf(data = hot_merged_2017, aes(fill = gly_per_land_low), linetype=0) +
  ggtitle("Glyphosate per Land Area", 
          subtitle="Hot Counties") +  theme_void(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(colour="gray25", size=14,hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks = element_blank(),
    rect = element_blank()
  ) +
  scale_fill_gradientn("Glyphosate usage",
                       colors = met.brewer("Tam"),
                       na.value = "grey50")
plot(hot_plot)

hot_leaflet <- hot_merged_2017 %>%
    st_transform(crs = "+init=epsg:4326") %>%
    leaflet(width = "100%") %>%
    addProviderTiles(provider = "CartoDB.Voyager") %>%
    addPolygons(
             stroke = FALSE,
             smoothFactor = 0,
             fillOpacity = 0.5,
             fillColor= ~
               colorNumeric("YlOrRd",gly_per_land_low)(gly_per_land_low)) %>%
  addLegend(pal = colorNumeric("YlOrRd", hot_merged_2017$gly_per_land_low),  # Define color palette
        values = hot_merged_2017$gly_per_land_low,  # Values to map
        title = "Glyphosate Usage (kg/km^2)",  # Legend title
        position = "bottomright"  # Position of the legend
    )
hot_leaflet

Change in Hot Counties Over Time

hot_merged_all <- list(hot_merged_2013, hot_merged_2014, hot_merged_2015,
                              hot_merged_2016, hot_merged_2017) %>% 
  rbindlist(fill=TRUE) %>% mutate(fips=as.numeric(fips))

hot_merged_sums <- hot_merged_all %>% group_by(YEAR) %>% 
  summarise(states=length(unique(STATEFP)),
            counties=length(unique(fips)),
            tmean85_fullyear = mean(tmean85_fullyear),
         tmean85_spray = mean(tmean85_spray),
         gly_low = sum(EPEST_LOW_KG, na.rm=TRUE),
         gly_high = sum(EPEST_HIGH_KG)) %>%
  filter(!is.na(YEAR))

hot_gly_use <- ggplot(hot_merged_sums) +
  aes(x = YEAR, y = gly_low) +
  geom_area() +
  geom_line(aes(x = YEAR, y = counties), colour = "#461711") +
  labs(
    x = "Year",
    y = "Glyphosate usage (kg)",
    title = "Glyphosate Usage 2013-2017",
    subtitle = "Hot Counties"
  ) +
  theme_bw() +
  theme(
    plot.title = element_text(hjust = 0.5),
    plot.subtitle = element_text(hjust = 0.5)
  )

plot(hot_gly_use)

Characterizing population in these areas

#put together simple summary table that counts how many people are in these hot counties. Include SVI & subscores.
hot_merged_2013 %>% group_by(STATEFP)  %>%

# I don't actually have data on population!

#Shiny App Good resources for helping learn Shiny: https://mastering-shiny.org/scaling-functions.html?q=map#functional-programming and example code creating an interactive map using Shiny: https://github.com/fverkroost/RStudio-Blogs/blob/master/interactive_worldmap_shiny_app.R These could also be useful: https://forum.posit.co/t/help-creating-a-simple-leaflet-interactive-plot-in-r-shiny/182374/2 https://rviews.rstudio.com/2019/10/09/building-interactive-world-maps-in-shiny/

# could create a binary data structure where there are 5 columns for each year and 448 rows for the compounds
#useful functions are compounddata_formap() and hot_compounddata_formap
ui <- fluidPage(
    # Application title
    titlePanel("Pesticides and Heat (2013-2017)"),
    
    fluidRow(
    column(3, selectizeInput("pesticide", label = "Pesticide", choices = usgs_compounds)), 
    column(3, selectInput("year", label = "Year", choices = c(2013, 2014, 2015, 2016, 2017))),
    column(4, numericInput('temp', label="85th %ile temp (C) in spray season", value=25)),
    column(3, checkboxInput('hot_counties', label= "Hottest counties only"))),
    
    # Application layout
    leafletOutput("map")
)

server = function(input, output, session) {
 # render the map
    output$map <- renderLeaflet({
      if(input$hot_counties==TRUE) {
        hot_leaflet_low(input$pesticide, input$year, input$temp)
      }
      else {
        all_leaflet_low(input$pesticide, input$year)
      }
    })
}

shinyApp(ui, server)
# look at Ag Census and ag worker density
# flexdash flexdashboard
# are there systematic patterns of more being applied when it is hot?
#make scatterplots of heat x pesticide where each location is an observation
saveRDS(usgs_ests, file="usgs_ests.rds")
saveRDS(usgs_compounds, file="usgs_compounds.rds")
saveRDS(us_county_2013, file="us_county_data.rds")
saveRDS(prism_county_2013_85s, file="prism_county_2013_85s.rds")
saveRDS(prism_county_2014_85s, file="prism_county_2014_85s.rds")
saveRDS(prism_county_2015_85s, file="prism_county_2015_85s.rds")
saveRDS(prism_county_2016_85s, file="prism_county_2016_85s.rds")
saveRDS(prism_county_2017_85s, file="prism_county_2017_85s.rds")
saveRDS(svi_2014, file="svi_2014.rds")
saveRDS(svi_2016, file="svi_2016.rds")
saveRDS(svi_2018, file="svi_2018.rds")

Functions

#-----functions----

# Show the names of all functions defined in the .Rmd
# (e.g. loaded in the environment)
lsf.str()
## all_leaflet_high : function (compound, year)  
## all_leaflet_low : function (compound, year)  
## compounddata_formap : function (compound, year)  
## hot_compounddata_formap : function (compound, year, temp_c)  
## hot_leaflet_high : function (compound, year, temp_c)  
## hot_leaflet_low : function (compound, year, temp_c)  
## map_svi : function (data, svi_var, year, percentile)  
## plot_gly : function (year, est_type)
# Show the definitions of all functions loaded into the current environment  
lsf.str() %>% set_names() %>% map(get, .GlobalEnv)
## $all_leaflet_high
## <srcref: file "" chars 26:21 to 49:1>
## 
## $all_leaflet_low
## <srcref: file "" chars 1:20 to 24:1>
## 
## $compounddata_formap
## <srcref: file "" chars 1:24 to 48:1>
## 
## $hot_compounddata_formap
## <srcref: file "" chars 50:28 to 52:1>
## 
## $hot_leaflet_high
## <srcref: file "" chars 76:21 to 99:1>
## 
## $hot_leaflet_low
## <srcref: file "" chars 51:20 to 74:1>
## 
## $map_svi
## <srcref: file "" chars 2:12 to 32:1>
## <bytecode: 0x55b3408941a8>
## 
## $plot_gly
## <srcref: file "" chars 1:13 to 52:1>
## <bytecode: 0x55b342a96940>

Versions

Information about the version of R and the packages using are below:

sessionInfo()
## R version 4.3.3 (2024-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.6 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] shiny_1.9.1        leaflet_2.2.2      terra_1.7-78       data.table_1.16.2 
##  [5] sf_1.0-18          here_1.0.1         ggpubr_0.6.0       MetBrewer_0.2.0   
##  [9] lubridate_1.9.3    forcats_1.0.0      stringr_1.5.1      dplyr_1.1.4       
## [13] purrr_1.0.2        readr_2.1.5        tibble_3.2.1       ggplot2_3.5.1     
## [17] tidyverse_2.0.0    tidyr_1.3.1        RColorBrewer_1.1-3
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.5            xfun_0.48               bslib_0.8.0            
##  [4] htmlwidgets_1.6.4       rstatix_0.7.2           tzdb_0.4.0             
##  [7] leaflet.providers_2.0.0 crosstalk_1.2.1         vctrs_0.6.5            
## [10] tools_4.3.3             generics_0.1.3          proxy_0.4-27           
## [13] fansi_1.0.6             highr_0.11              pacman_0.5.1           
## [16] pkgconfig_2.0.3         KernSmooth_2.23-22      lifecycle_1.0.4        
## [19] farver_2.1.2            compiler_4.3.3          munsell_0.5.1          
## [22] codetools_0.2-20        carData_3.0-5           httpuv_1.6.15          
## [25] htmltools_0.5.8.1       class_7.3-22            sass_0.4.9             
## [28] yaml_2.3.10             Formula_1.2-5           later_1.3.2            
## [31] pillar_1.9.0            car_3.1-3               jquerylib_0.1.4        
## [34] classInt_0.4-10         cachem_1.1.0            abind_1.4-8            
## [37] mime_0.12               tidyselect_1.2.1        digest_0.6.37          
## [40] stringi_1.8.4           labeling_0.4.3          rprojroot_2.0.4        
## [43] fastmap_1.2.0           grid_4.3.3              colorspace_2.1-1       
## [46] cli_3.6.3               magrittr_2.0.3          utf8_1.2.4             
## [49] broom_1.0.7             e1071_1.7-16            withr_3.0.1            
## [52] promises_1.3.0          scales_1.3.0            backports_1.5.0        
## [55] timechange_0.3.0        rmarkdown_2.28          ggsignif_0.6.4         
## [58] hms_1.1.3               evaluate_1.0.1          knitr_1.48             
## [61] viridisLite_0.4.2       rlang_1.1.4             Rcpp_1.0.13            
## [64] xtable_1.8-4            glue_1.8.0              DBI_1.2.3              
## [67] rstudioapi_0.16.0       jsonlite_1.8.9          R6_2.5.1               
## [70] units_0.8-5